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Abstract 

The  aim  of  this  report  is  to  present  a  consistent  theory  for  the  deformation  of  a 
naturally  curved  and  twisted  anisotropic  beam.  The  proposed  formulation  natu¬ 
rally  extends  the  classical  Saint- Venant  approach  to  the  case  of  curved  and  twisted 
anisotropic  beams.  The  mathematical  model  developed  under  the  Eissumption  of 
span- wise  uniform  cross-section,  curvature  and  twist,  can  take  into  account  any 
kind  of  elastic  coupling  due  to  the  material  properties  and  the  curved  geometry. 
The  consistency  of  the  math-model  presented  and  its  generality  about  the  cross- 
sectional  shape,  make  it  a  useful  tool  even  in  a  preliminary  design  optimization 
context  such  as,  for  exeunple,  the  aeroelastic  tailoring  of  helicopter  rotor  blades. 

The  advantage  of  the  present  procedure  is  that  it  only  requires  a  two-dimensionail 
discretization;  thus,  very  detailed  malyses  can  be  performed  and  interlamineir 
stresses  between  laminae  can  be  evaluated.  Such  analyses  would  be  extremely  time 
consuming  if  performed  with  stamdard  finite  element  codes:  that  prevents  their 
recursive  use  as  for  example  when  optimizing  a  beam  design. 

Moreover,  as  a  byproduct  of  the  proposed  formulation,  one  obtains  the  consti¬ 
tutive  law  of  the  cross-section  in  terms  of  stress  resultemt  and  moment  and  their 
conjugate  strain  measures.  This  constitutive  law  takes  into  account  any  kind  of  elas¬ 
tic  couplings,  e.g.  torsion- tension,  tension-shear,  bending-shear,  and  constitutes  a 
fundamental  input  in  aeroelastic  suialyses  of  helicopter  blades. 

Four  simple  examples  are  given  in  order  to  show  the  principal  features  of  the 
method. 
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1  statement  and  Scope 

The  relevance  of  aeroelastic  tailoring  of  helicopter  blades  is  very  well  known  in  rotor- 
craft  community.  In  fact,  the  elastic  couplings  can  strongly  influence  the  aeroelastic 
behavior  of  rotor  systems  and  even  control  their  stability.  Stability  analyses  of  heli¬ 
copter  rotors  axe  often  performed  using  speciad  beam  models  to  simulate  the  blade 
dynamic  behavior,  i.e.  the  blade  is  considered  ais  a  one-dimensional  continuum 
with  a  general  form  of  elaistic  couplings.  Therefore,  the  constitutive  law  expressed 
in  terms  of  stress  resultant  and  moment  on  the  cross-section  and  their  conjugate 
strain  meaisures  is  a  fundamental  input  for  aeroelastic  anaJyses.  In  this  regaird,  the 
elaistic  couplings,  such  as  tension- torsion,  tension-shear,  bending-sheax,  etc.,  axe  of 
fundamentad  importance  in  performing  aeroelaistic  tauloring. 

Such  couplings  arise  from  two  different  sources,  namely  materiad  properties  and 
geometric  shape  of  the  blade.  Obviously,  a  very  generad  way  in  order  to  compute 
the  elastic  constitutive  relations  is  using  a  general  three-dimensional  finite  element 
code,  but  a  drawback  of  this  approach  is  the  amount  of  computer  time  required 
that  cannot  be  accepted  during  the  preliminairy  design  phases  in  which  different 
fiber  orientations  and  materials  axe  to  be  analyzed  in  order  to  meet  the  design 
requirements. 

An  adternative  approach,  in  order  to  circumvent  this  problem,  is  to  take  some 
suitable  simplifying  assumptions  allowed  by  the  blade  geometry  aind  by  the  aims 
of  the  preliminairy  design  itself,  as  for  example  the  constancy  of  the  cross-section 
properties  along  the  beam  axis.  This  assumption  leads  to  a  dramatic  improvement 
of  the  analysis  since  the  three-dimensional  problem  can  be  eaisily  reduced  to  a  two- 
dimensional  one.  In  fact  the  object  of  the  ainalysis  is  now  reduced  to  a  simple  beam 
slice  instead  of  the  entire  blade.  By  the  use  of  this  simplified  physical-mathematical 
model,  the  eleistic  properties  of  the  blade  cross-section  can  be  optimized  emd  prop¬ 
erly  tailored  and  only  few  three-dimensional  analyses  are  required  to  eissess  the  final 
design. 

Solutions  in  closed  form  of  composite  cross-sections  are  not  available  and  a 
two-dimensional  finite  element  model  of  the  blade  section  is  required.  This  finite 
element  model  is  a  special  one  since  it  is  able  to  recover  the  three-dimensional 
solution  with  the  only  assumption  of  constancy  of  the  cross-section.  This  kind 
of  einalysis,  developed  by  the  writer  and  its  colleagues,  was  already  available  for 
straight  beeuns.  The  sum  of  this  work  is  to  extend  this  approach  to  naturally  twisted 
and  curved  beams  retaining  the  assumption  of  the  constancy  of  the  cross-section 
properties.  The  influence  of  the  pre-twist  become  mandatory  in  the  analysis  of  tilt 
helicopter  rotors  like  the  JVX,  in  which  the  amount  of  twist  required  in  airplane 
mode  by  aerodynamic  consideration  is  much  larger  than  that  usually  required  for 
conventional  helicopter  rotors:  therefore  a  stronger  influence  can  be  eJso  expected 
from  a  structural  point  of  view. 
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The  research  work  reported  herein  can  be  summarized  as  follows; 

•  Development  of  the  mathematical  formulation  of  the  mechanics  of  a  beam 
slice,  b2ised  on  the  principle  of  virtual  work  amd  taking  into  account  the  ef¬ 
fects  of  a  constant  pre- twist  and  pre-bending.  The  possibility  of  a  cross-section 
tilted  with  respect  to  the  axis  has  been  investigated.  This  section  extends  the 
approach  already  available  for  straight  beams.  One  result  of  such  an  anadysis  is 
the  solution  of  the  cross-section  problem  in  terms  of  displacement  and  stress 
distribution  due  to  six  independent  and  equilibrated  load  conditions  corre¬ 
sponding  to  the  components  of  the  cross-section  stress  resultamt  and  moment. 
Moreover,  as  a  byproduct  of  this  analysis,  the  six  by  six  compliance  matrix 
of  the  cross-section  is  obtained.  In  addition,  the  eigenvalue  analysis  in  terms 
of  warping  displacements  is  performed  and  the  extremity  solutions  computed. 
The  results  of  these  analyses  can  be  used  as  input  in  subsequent  dynamic 
analyses  in  which  the  beam  is  treated  zis  a  one-dimensional  continuum. 

•  Formulation  of  different  finite  elements  in  anisotropic  materiad  like  the  isopara¬ 
metric  plane  element,  the  lamina  element,  the  isoparametric  panel  and  the 
stringer  element. 

•  Modification  of  the  pre-existing  computer  code  for  straight  beam,  in  order  to 
include  the  curved  amd  tilted  geometry.  Description  and  organization  of  the 
computer  code  based  on  the  finite  element  method. 

•  Development  of  some  numerical  test  caises  and  comparison  with  NASTRAN 
results. 
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Theoretical  Development 

2.1  Introduction 

The  design  of  composite  bl2ides  of  helicopter  rotors  demands  the  analysis  of  three- 
dimensionzkl  stress  states  including  interlaminar  stresses.  Despite  the  power  of  mod¬ 
ern  computers,  standeird  three-dimensional  finite  element  approximations  of  the  en¬ 
tire  rotor  blade  are  not  yet  considered  feasible,  because  of  the  huge  computer  effort 
required  to  achieve  a  reasonable  degree  of  accuracy  in  modeling  the  material  proper¬ 
ties  of  the  blade  cross-section  (Hodges,  1990b).  In  fact  when  dealing  with  composite 
beams,  particular  care  must  be  taken  in  order  to  model  the  material  properties,  es¬ 
pecially  when  interlaminar  stresses  must  be  evaluated.  The  example  of  the  circular 
composite  tube,  see  Fig.  1  (from  Giavotto  et  a/.,  1983),  shows  that  the  stress  state 
is  three-dimensional  even  under  pure  torsion,  hence  a  detailed  model  of  the  cross- 
section  is  essential  in  order  to  correctly  capture  it.  Moreover,  the  stacking  sequence 
of  different  laminae  may  have  a  strong  infiuence  on  the  cross-section  distortion  and 
on  the  stress  distribution,  as  shown  for  example  by  Ghiringhelli  and  Sala  ( 1990);  as 
a  simple  example  refer  to  Fig.  2  (from  Ghiringhelli  and  Sala,  1986),  where  the  in¬ 
plane  section  distortions  of  a  flat  specimen  loaded  in  uniform  tension  cire  shown  for 
two  different  cross-ply  laminates.  These  examples  stress  that  for  composite  beams 
not  only  is  the  out-of-plzme  distortion,  usually  accounted  for,  significant,  but  adso 
the  in-plane  distortion  must  be  taken  into  account  in  order  to  correctly  model  the 
stress  behavior.  It  is  importzuit  to  observe  as  well,  that  the  interlaminar  stresses 
between  contiguous  laminae  must  be  continuous  and,  since  the  material  properties 
are  in  general  different,  the  conjugate  strains  are  discontinuous.  This  fact  consti¬ 
tutes  a  prerequisite  of  any  displacement-based  finite  element  approximation  and,  in 
order  to  model  this  strain  discontinuity,  the  element  size  must  not  be  greater  thtui 
the  thickness  of  the  lamina. 

The  high  degree  of  detail  required  to  model  such  problems  cannot  be  tolerated 
in  a  direct  three-dimensional  approach,  since  the  high  number  of  degrees  of  freedom 
prevents  its  use  in  preliminwy  design  pheises  and  in  am  optimization  context.  On 
the  contrary,  it  would  be  better  to  face  the  problem  in  two  consecutive  steps  deal¬ 
ing  separately  with  the  cross-sectional  analysis  and  with  the  global-beaun  analvsis. 
In  the  first  step  the  three-dimensional  stress  analysis  of  the  beam  cross-se-tion, 
modeled  on  a  two-dimensional  domain,  is  performed  leading  to  a  differential  prob¬ 
lem  with  respect  to  the  span-wise  coordinate  of  the  beam.  This  step  is  generally 
performed  with  a  finite  element  procedure,  as  proposed  by  Giavotto  et  al.  (1983), 
although  an  analytical  approach  giving  solutions  in  closed  form  could  fit  as  well  for 
very  simple  cross-sections.  In  this  step  particular  solutions  under  prescribed  stress 
resultants  are  obtained,  giving  the  stiffness  of  the  cross-section  and  its  generalized 
warping,  i.  e.  section  out-of-plane  and  in-plane  distortion.  Moreover  eigensolutions 
can  be  obtained  giving  the  diffusion  length  of  self-equilibrated  modes,  that  can  be 
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Figure  1:  Thin-walled  composite  cylinder  [-I-45/-45]  loaded  in  torsion  (not  to  scale). 

superimposed  to  the  central  solutions  to  account  for  extremity  effects,  either  at  the 
beam  ends  or  in  the  neighborood  of  concentrated  loads.  These  modes,  zis  p>ointed 
out  by  Rehfield  et  al.  (1990),  could  also  be  used  eis  additional  kinematical  variables 
when  modeling  thin-walled  beams,  while  their  role  in  axially  compressed  beams  and 
pzmels  hsis  been  dealt  with  by  Merlini  (1988). 

The  second  step  is  mainly  devoted  to  the  behavior  of  the  entire  beam,  this  in 
generjd  being  naturally  curved  and  twisted:  here  the  beEim  is  considered  as  a  one- 
dimensioned  continuum  zmd  its  constitutive  law  is  taken  from  the  previous  step. 
This  kind  of  approach  naturally  extends  the  well  known  Saint- Venant  approach  to 
pretwisted  zmd  curved  anisotropic  beams,  and  it  is  gmning  more  and  more  attention 
in  rotor  blade  modeling,  as  discussed  in  a  review  by  Hodges  (1990b).  Moreover,  this 
approach  seems  to  be  suitable  for  geometrically  nonliner  beam  juialyses,  both  static 
and  dynamic,  where  strains  can  be  assumed  to  be  small  amd  nonlinearities  confined 
within  the  one-dimensional  beam  model,  as  shown  by  Borri  and  Mamtegazza  (1985) 
and  Borri  and  Merlini  (1986).  This  idea  can  be  found  in  the  papers  by  Parker 
(1979a,  b),  who  followed  an  asymptotic  procedure  and  used  the  Saint- Venamt  wairp- 
ing  function,  and  by  Berdichevsky  (1981)  and  Berdiclievsky  and  Stairoselsky  (1983), 
who  proposed  a  vau'iationad-asymptoticaJ  method.  A  generad  aissessment  of  this  ap¬ 
proach  hats  been  set  up  by  Hodges  amd  his  co-workers  and  is  plainly  described  in 
some  recent  papers:  Hodges  (1990a),  Atilgan  and  Hodges  (1990),  Atilgan  ti  al. 
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Figure  2:  Flat  composite  specimen  loaded  in  social  tension  and  in-plane  distortion 
for  two  different  stsw:king  sequencies. 
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(1991)  and  Hodges  ei  al.  (1991a).  An  application  to  beam  free- vibration  analysis 
can  be  found  in  Hodges  et  al.  (1991b). 

It  is  not  in  the  purpose  of  this  report  to  review  the  intensive  work  towards 
comprehensive  beam  modeling  found  in  the  literature:  an  excellent  review  is  found 
in  Hodges  (1990b).  Nevertheless,  at  least  the  following  papers  must  be  mentioned 
for  their  significant  contribution  to  the  one-dimensional  formulation  and  analysis  of 
the  nonlinear  behavior  of  either  rectilinear  or  space-curved  beams:  Reissner  (1973, 
1981),  Simo  (1985),  Simo  and  Vu-Quoc  (1988),  Hodges  (1987a,  b),  Danielson  and 
Hodges  (1987,  1988),  Bauchau  and  Hong  (1987),  Cardona  and  Geradin  (1988)  and 
lura  and  Atluri  (1988,  1989).  Other  pai>ers  are  mentioned  here  contributing  to  the 
three-dimensional  analysis  of  the  beam  cross-section.  Theories  by  Bauchau  (1985) 
2ind  Bauchau  zind  Hong  (1988)  are  demonstrated  to  behave  well  but  the  assumption 
of  indeformability  of  the  cross-section  in  its  own  plzme  seems  too  strong  for  a  general- 
purpose  zuialysis  of  comp>osite  beams.  Kosmatka  and  Dong  (1991)  proposed  an 
analytical  model  of  the  anisotropic  cross-section  yielding  the  global  properties  of 
the  section,  but  their  work  is  restricted  to  homogeneous  beams.  Finally,  the  work 
by  Stemple  and  Lee  (1988, 1989)  is  mentioned,  though  it  is  not  relevzint  to  einy  of  the 
two  steps  stated  above:  they  developed  a  peu'ticulair  fully  three-dimensional  finite 
element  procedure  accounting  explicitly  for  warping  of  anisotropic  sections,  but  it 
seems  that  the  prerequisite  of  stress  continuity  and  strain  discontinuity  between 
laminae  could  hardly  be  satisfied  in  practice  by  their  approach. 

The  present  report  aims  to  contribute  to  the  first  step  outlined  above.  A  quite 
general  theory  is  formulated  for  the  three-dimensional  cross-section  analysis  ac¬ 
counting  for  initial  twist  and  curvature  of  anisotropic  and  nonhomogeneous  beams. 
The  constitutive  law  of  the  section  obtained  this  way  could  then  be  used  in  most 
of  the  one-dimensional  theories  for  space-curved  beams  proposed  by  many  authors, 
such  as  those  mentioned  above.  The  generality  of  the  proposed  theory  is  only  lack¬ 
ing  in  the  sense  that  it  allows  for  constant  twist  eind  curvature  along  the  be^un 
span;  from  a  practical  point  of  view,  this  restriction  can  be  suitably  overcome  in 
the  one-dimensionzd  analysis  by  means  of  a  spatial  interpolation,  using  beam  finite 
elements.  Although  this  theory  was  motivated  by  the  demands  inherent  to  mod¬ 
eling  helicopter  rotor  blades,  it  is  believed  that  it  constitutes  a  valid  tool  in  many 
other  structured  fields.  Thus  it  would  be  proper  to  begin  to  state  what  a  beam  is. 
It  is  worth  noting  that  what  is  usuedly  called  a  beam  is  not  so  univocally  defined 
since  there  exist  many  different  approximations.  Therefore  let  us  recedl  some  beisic 
definitions  emd  restrictions. 

There  are  at  least  two  fundamented  restrictions  for  a  structured  member  to  state 
that  it  is  a  beeun:  one  involves  only  geometriced  considerations,  while  the  other 
concerns  structured  mechanics.  From  a  geometrical  point  of  view  a  beam  cem  be 
defined  ets  a  solid  treK:ed  by  the  rigid  motion  of  a  cross-section.  One  particuleu:  point 
of  this  section  can  be  taken  as  the  reference  point  and  the  line  traced  by  this  point 
as  reference  line,  which  is  required  to  possess  a  certain  level  of  continuity.  The 
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geometric  beam  concept  also  includes  some  slenderness  restrictions  mainly  related 
to  the  width  of  the  cross-section  that  must  be  much  smaller  than  the  length  of  the 
reference  line.  Moreover,  if  the  cross-section  is  not  constant,  its  vziriation  should  be 
restricted  to  be  only  moderate. 

Prom  the  structural  point  of  view,  in  addition  to  the  influence  of  the  boundary 
conditions,  we  must  mention  at  least  two  fundamental  assumptions,  one  being  rele¬ 
vant  to  the  material  properties  and  the  other  to  the  applied  loads:  these  hypotheses 
concern  the  eixial  distribution  of  the  material  properties  and  the  applied  loads  that 
must  be  very  smooth.  As  a  consequence,  every  structural  variable  like  displacement 
or  stress  should  have  a  gradient  along  the  beam  zocis  that  is  much  lower  than  the 
gradient  along  the  cross-section  coordinates. 

Now,  it  should  be  understood  that  if  a  solid  must  behave  like  a  beam,  the 
restrictions  outlined  above  must  hold  not  only  before  the  deformation  but  also 
during  the  deformation  and  in  the  deformed  configuration  too.  This  issue  leads  to 
the  concept  of  a  reference  cross-section  which  could  be  tilted  with  respect  to  the 
reference  line,  and  even  not  exactly  planar,  in  the  undeformed  configuration  as  well. 
Furthermore,  this  concept  adlows  for  an  easy  implementation  of  updated- Lagrangean 
incremental  techniques.  All  these  hypotheses  allow  us  to  consider  a  beam  as  a  one¬ 
dimensional  continuum,  a  point  of  which  reacts  against  the  deformation  with  one 
resultant  and  one  moment,  about  the  reference  point,  of  the  stresses  acting  upon 
the  cross-section.  With  the  previous  restrictions  in  mind,  these  stress  resultant  and 
moment,  for  a  specified  reference  line,  should  be  independent  of  the  orientation  of 
the  reference  cross-section,  which  obviously  cannot  be  too  far  from  the  normal  one. 
Moreover,  some  freedom  should  also  hold  for  the  choice  of  the  reference  line:  in 
fact,  if  we  think  about  a  beam  which  is  naturally  twisted,  a  section  will  be  normal 
to  only  one  reference  line  and  it  will  not  be  normal  to  any  other  reference  line. 

Even  if,  in  our  opinion,  the  main  difficulties  related  to  this  subject  are  not 
confined  to  the  development  of  a  consistent  mathematical  model,  but  instead  in  the 
correct  identification  of  the  elastic  properties  of  composite  materials,  it  still  remeiins 
difficult  to  give  a  comprehensive  definition  of  a  beam  and  we  can  accept  many  of 
them,  depending  on  a  good  judgment  of  what  the  goal  is  and  what  we  can  use  to 
obtain  it. 

2.2  One-dimensional  Beam  Geometry 

In  a  prescribed  fixed  orthogonal  reference  frame  (O,  £?,)  (i  =  1,2,3),  let  x  be  the 
position  of  a  generic  point  P  of  a  reference  line  I  and  let  t  be  the  tangent  to  I  at  P. 
We  cam  write 

x  =  x(s)  =  (1) 


where  s  denotes  the  mc  length  along  /. 


For  simplicity  let  A{s)  be  a  plane  cross- section.  Associated  to  it  define  sm 
orthogonal  triad  (P,  ej:  the  unit  vectors  ei(s)  and  62(5)  belong  to  the  cross-section 
while  63(5)  =  ei(s)  x  62(5)  may  be  tilted  with  respect  to  t{s);  however,  we  assume 
that  €3(3)  •  t{s)  >  0  holds  true  anywhere.  The  orientation  of  the  cross-section  can 
be  specified  through  a  rotation  of  the  section  at  s  =  0, 


Ci(s)  =  fi(r)-e,(0)  (2) 

where  r(s)  is  the  rotation  vector  and  jR(r)  is  the  corresponding  rotation  tensor 
which  by  definition  satisfies  the  orthonormality  property  R{r)  •  R^{r)  =  I.  It  is 
well  known  that  the  last  can  be  given  the  following  form 


R(r )  =  /  -f  sin  (^iit  X  /  -f  ( 1  —  cos  (i>)k  X  k  X  I  (3) 

where  4>(s)  =  \Jr(s)  ■  r(s)  is  the  magnitude  of  the  rotation  vector  r{s)  =  (f>k  eind  k 
the  unit  vector  of  the  axis  of  the  rotation. 

The  generalized  curvature  vector  c(s)  of  the  beam  can  be  introduced  through 
the  derivative  of  the  rotation  that  can  be  cast  in  the  form 


dR{T) 

ds 


=  c(s)  X  R(r) 


due  to  the  orthogonality  property  of  the  rotation  itself.  From  eqn  2  it  follows  that 


de,(  s) 

~^  =  c(s)xc,(s).  (4) 

In  this  report  we  confine  ourselves  to  deal  with  the  special  case  of  a  helicoidal 
beam  for  which  the  following  limitations  hold: 

^(t(s)  •  e,(s))  =  0 

From  these  equations,  taking  into  account  eqn  4,  we  obtain 


dt(s) 

ds 


=  c{s)  X  t{s} 


dc{s) 

ds 


=  0 


so  that  the  curvature  c  is  constant  while  the  vector  t{s)  rotates,  i.e.  it  can  be 
expressed  in  the  following  form: 


t(s)  =  i2(r)  •  t(0). 


(5) 


off*  f  3  ) 

Moreover,  since  in  this  case  it  can  be  shown  that  c  =  — it  follows  that 

ds 


r{s)  =  sc  =  sck 


<f>(s)  =  sc. 


(6) 
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i.e.  the  axis  of  rotation  is  constant  and  coincides  with  the  axis  of  the  curvature. 

For  such  a  helicoidal  beam  the  equation  of  the  line  /  can  be  written  as  the 
following: 

aj(s)  -  a  =  <^(s)afc  +  H(r)  •  (®(0)  -  a)  (7) 

c  c 

and  A  is  an  indeterminate  scalar  quantity  which  can  arbitrarily  be  set  to  zero  by 
choosing  the  vector  a  such  that  *(0)— a  be  orthogonal  to  Jb.  Equation  7  corresponds 
to  the  so  called  helicoidal  decomposition  of  a  rigid  displacement  into  a  rotation  zmd 
a  translation  parallel  to  it.  The  tangent  vector  t(0)  cam  be  resolved  in  the  following 
way, 

t(0)  =  fc  <8i  fc  •  t(0)  +  il-k®k)- 1(0)  = 

=  c{ak  +  k  X  {x{Q)  -  a))  (8) 

where  the  identity  k  X  k  X  I  =  k®k  —  I  has  been  taken  into  account.  From  eqn  5, 
taking  eqns  8  and  7  into  account,  we  obtain  the  expression  for  the  tangent  vector 
t{s)  to  the  reference  line  /: 

t{s)  =  c{ak  +  k  X  (a;(s)  —  a)) . 

Equation  7  can  be  easily  obtained  by  integrating  eqn  1  with  eqn  5  on  the  curvi¬ 
linear  abscissa  s.  At  first  the  tensor  S(r)  is  introduced, 


S(r)  =  -  r  R{r)ds  =  1  + 
s  Jo 

leading  to 


1  -  cos  ,  , ,  sin  ,  ,  _ 

k  X  I  +  {1 - —)k  X  k  X  I 


x{s)  =  x{0)  -f-  sS{r)  •  t(0). 

Then  recalling  eqns  8  and  6  and  taking  the  following  property  into  account, 

Rir)  =  I  +  S{r)  r  X  I 


eqn  7  is  obtained. 

It  is  worthwhile  noting  that  this  description  of  helicoidal  beams  is  invariant 
with  respect  to  the  choice  of  the  reference  line.  In  order  to  clear  up  this  concept, 
let  /*  denote  a  new  reference  line  different  from  I,  and  let  x‘{s)  and  x{s)  be  the 
position  vectors  of  two  points  P*  and  P  belonging  to  those  lines  and  to  the  same 
cross-section.  Hence  the  line  /*  can  be  defined  through  the  equation 

a!*(s)  =  *(s) -f- 60(3)77“  (9) 

with  77“  (a  =  1,2)  constant.  Recalling  eqn  2  we  cam  write 

x*(s)  -  x(3)  =  H(r)  •  (x’(0)  -  x(0)). 
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Then  eqn  7  transforms  to 

x’(s)  —  a  =  4>{s)ak  +  R{r)  ■  (x*(0)  —  a). 

In  order  to  completely  change  the  reference  line  we  should  also  change  the  indepen¬ 
dent  variable  to  the  new  arc  length  s*  on  /*.  By  definition,  we  have 


ds*  =  ^dx‘{s)  •  dx‘{s)  =  jds 

and  since  in  the  pr^ent  case  the  the  Jacobian  j  is  constant,  implying  s‘  =  js, 
from  eqn  6  we  can  obtain  the  new  curvature  c*  =  ^.  Hence,  by  the  definition 


t\s)  = 


dx*{s)  Idx’(s) 


,  we  easily  obtain  the  tangent  vector  to  the  new  reference 


ds"  j  ds 

line  /*: 

t’{s)  =  c'  {ak  +  k  X  (x’{s)  —  a)) . 

Although  the  tangent  and  the  curvature  vectors  depend  on  the  choice  of  the 
reference  line,  the  parameters  a  and  a  describing  the  elicoidal  beam  are  invariant: 
with  reference  to  the  new  line  they  can  be  written  as 

t’(0) 


a  =  k  ■ 


a  =  **(0)  +  k  X  ^  yf. 


c*  c’ 

where  A*  =  A  —  fc  •  (*‘(0)  —  x(0))  is  a  new  indeterminate  scalar. 


2.3  One-dimensional  Beam  Equations 


In  a  generic  deformed  position  of  the  beam,  let  T{s)  and  Af  (s)  denote  respectively 
the  resultant  and  the  moment  of  the  stresses  acting  upon  the  cross-section  .4(s). 
In  the  small- displacement  linear  theory,  the  equilibrium  equations  of  the  beam  are 
written  as 


dT(s) 

ds 


-  0 


dM{s) 

ds 


—  T(s)  X  t(s)  =  0 


(10) 


where  the  explicit  influence  of  the  distributed  external  forces  is  omitted,  since  our 
interest  is  confined  to  homogeneous  solutions.  These  equations  can  be  immediately 
integrated  leading  to: 


T{3)  =  T(0) 

M(s)  =  M(0)  +  T(0)  X  (x(s)  -  x(0)). 


(11) 


Both  of  eqn  11  have  a  very  well-known  and  clear  physical  meaning.  Looking  at  the 
expression  for  x(s)  —  *(0)  in  eqn  7  we  can  see  that  the  generaJ  solution  depends  on 
s  by  means  of  linear  and  circular  functions. 
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Moreover  it  is  interesting  to  obtain  the  solutions  of  eqn  10  in  terms  of  cross- 
sectional  components,  i.e.  Ti{s)  =  T{s)  ■  e,  and  Mi{s)  =  Mis)  •  Ci-  For  simplicity 
we  assume  the  fixed  reference  triad  to  be  oriented  as  the  cross-section  triad  at 
s  =  0,  so  that  E,  =  ei(0),  and  define  T(s)  =  Ti{s)Ei,  Mis)  =  i.e. 

T’(s)  =  •  T(s),  Mis)  =  il^(r)  •  Mis).  This  is  equivalent  to  pulling  back  the 

stress  resultant  and  moment  to  the  cross-section  at  s  =  0.  By  observing  that 

<^Tis)  ,  fdfis)  .  *  A 

where  c  =  iZ^(r)  ■  c  =  c,  the  one-dimensional  beam  equilibrium  equations  in  terms 
of  T(s)  and  Mis)  axe  easily  written 


dTjs) 

ds 

dMjs) 

ds 


+  cx  Tis)  =  0 
-f-  c  X  Mis)  +  ix  Tis)  =  0 


where  i  =  ir)  •  tis)  =  t(0).  These  equations  can  be  grouped  into  one  equation 
of  the  form 

dSis) 


ds 


+  T-Sis)  =  0 


(12) 


where 


Sis)  =  iTis).Mis)) 


r  = 


cx  I 

0 

i  X  I 

cx  I 

Observing  that  the  operator  T  is  constant  and  the  following  property  holds, 


-I-  2c't‘'  -F  =  0. 


as  it  can  be  easily  seen,  by  differentiating  eqn  12  we  obtain 


d^Sjs)  _  -4 
ds^ 


Sis)  and  — '  Sis),  whence 
as® 


cPS(s) 

ds2 


r 


Sis), 


d^Sjs)  .d^Sjs)  ,<PSis) 
ds^  ds'*  ^  ds^ 


=  0. 


(13) 


The  peculiarity  of  this  equation  is  that  each  component  of  the  vector  <S(s)  is  gov¬ 
erned  by  the  same  differential  equation.  Thus,  the  eigenvalues  are  immediately 
found  to  be  the  solutions  of  the  characteristic  equation 

A^(A  —  tc)^(A  -f  ic)^  =  0  i  =  \/^. 
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As  a  consequence,  eqn  12  admits  a  solution  of  the  following  form; 

5(s)  =  (Ao  +  +  (Ac  +  <^Bc)  cos  <i>  +  +(A,  +  <f>Bs)sm<i>.  (14) 

Imposing  the  initial  condition  at  s  =  0,  by  means  of  subsequent  derivatives  of  eqn  12 
at  s  =  0  we  can  obtmn: 


Ao  = 


Ac  = 


A,= 


(j  +  2c-^t'‘ +  -5(0) 

•  (21  +  ■  5(0) 

1/20-^"*  •  (5/+  3c-2t")  •  5(0) 


=  -c-^r-( 


Bc  = 
Bs  = 


r  •  ^1  +  2c-2T*  +  ■  5(0) 

-\/2c-^r^  •  +  c-2T^)  •  5(0) 

-i/2c-2r"  •  +  c-^r')  •  5(0) 


Then  we  can  express  eqn  14  in  the  following  form: 


where 


5(s)  =  4-  Y^(-c)-'‘M<t>)T  J  •  5(0) 

Vi  =  0 


V’2  =  2  —  2  cos  0  —  sin 
03  =  ^<f>  cos  (^  —  I  sin  o  +  2<?1> 

04  =  1  —  cos  0  —  ^(f>  sin  (j) 

ii^5  =  0  +  ^(^cos  4>-  Isin  <t>. 

It  is  intersting  to  derive  an  approximation  of  these  expressions  for  small  curva¬ 
ture.  This  is  easily  obtained  observing  that: 

lim(-c)-V*  = 

Then  5(s)  can  be  approximated  by 

5(^)  =  +  -.5(0).  (17) 

te=I 

This  is  a  consistent  approximation  since  for  small  curvature  the  differential  equi- 

d®5(s) 

librium  eqn  13  is  approximated  by  — 7-; —  =  0,  which  admits  eqn  17  as  solution. 

ds° 

2 

Moreover,  since  for  c  =  0  also  'T  =0,  for  this  particular  case  we  have 

Ao  =  ^(0)  ^Bq  =  —s'T  •  5(0)  Ac  =  Aj  =  Be  —  B t  ~  0 
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so  that  we  obtain  a  linear  solution  for  the  generalized  stress  resultant  in  the  form 

S{s)  =  (/  -  st) .  5(0). 


2.4  Three-dimensional  Beam  Deformation 

In  order  to  obtain  the  constitutive  equations  of  the  cross-section,  i.e.  the  cross- 
sectional  rigidity  matrix,  let  us  consider  the  beam  as  a  three-dimensional  continuum. 

The  position  vector  xq  of  a  generic  material  point  Q  in  the  cross-section  can  be 
expressed  as 

=  x(3) -I- 

where  x(s)  is  the  position  of  the  point  P  on  I  and  ^“(q  =  1,2)  are  the  czurte- 
sian  coordinates  along  the  cross-sectional  axes  eo,(s).  This  representation  induces 
curvilinear  coordinates  (A:  =  1,2, 3),  =  s,  with  covariant  basis  vectors 


=  Co  g3  =  t  +  c  X  b 


where  b  =  Xq  —  x  is  the  vector  from  point  P  to  point  Q.  Moreover,  denoting  by  g 
the  determinant  of  the  metric  tensor  such  that  y^g  =  Qi  ^  92  '  Osi  contravariant 
basis  vectors  are  given  by 


<73  X  C3  X  Ca 

g  - - - - 


9^  =  —- 

^  n/? 


In  the  following,  in  order  to  compute  the  strain  tensor,  we  need  the  expression  for 
the  displacement  gradient  tensor.  Let  Vs  denote  the  gradient  of  the  displacement 

,  ds 

vector  s,  i.e.  Vs  =  9  ^dS  denote  the  covetrieint  derivative  of  s  along 

the  direction  d,  i.e.  V^s  =  <i  •  Vs.  In  the  same  way  the  covariant  derivatives  of  s 
along  the  cross-sectional  basis  vectors  au’e  denoted  by  Vg.s  and  are  given  by 


Ve.s  =  e,  Vs  =  a  g* 


ds 


Substituting  in  this  equation  the  expressions  for  the  contravariant  basis  vectors  we 
easily  obtain 

Ve„s 
VcsS 

vy 

ds 

,  d°8  _  ds 

^“cre  —  c  X  8  denotes  the  corotationzil  derivative  with  the  cross- 

sectional  triad,  with  respect  to  the  eu-c  length. 


ds 

d^^ 

1  I 

(ds 

ds 

[de 

1  1 

(&>s 

[d^ 

-f  C  X  S  -  ea 

Within  the  approximation  of  small  displacements  the  straiin  tensor  is  the  sym¬ 
metric  part  of  the  displacement  gradient  tensor,  i.e.  e  =  |  (Vs  and  there 

is  no  need  to  distinguish  among  the  various  nonlinear  stress  tensors.  Then  let  <r 
denote  the  stress  tensor  at  the  generic  point  Q,  corresponding  to  the  displacement 
field  s  zmd  let  us  assume  the  validity  of  Hooke’s  law,  i.e.  er  =  D  •  €  where  D  is 
the  elasticity  tensor.  We  do  not  specify  any  restriction  to  the  nature  of  the  elastic 
tensor  D  except  that  it  is  symmetrical  and  jxjsitive  definite  and  its  corotational 
derivative  is  zero.  Moreover,  let  n  denote  the  unit  normal  to  the  cross-section  and 
p  =  •  n  the  corresponding  stress  vector. 

For  the  sake  of  simplicity,  we  assume  that  the  externJil  loads  are  applied  only 
at  the  end  sections  of  the  beam,  and  the  distributed  loads  are  negligible  or  absent. 
For  our  purposes,  these  hypotheses  do  not  affect  the  results  we  are  looking  for. 


2.5  Virtual  Work  Principle  and  F.E.M.  Approximation 

Using  the  notation  of  the  previous  paragraphs,  the  principle  of  virtual  work  written 
for  an  infinitesimal  beam  slice  reads 

J  tr{6€  ■  D  •  €)y/gdA  =  ^ 

where  6s  denotes  the  virtual  displacement  of  the  generic  point  Q  and  = 
I  (ySs  -f  is  the  corresponding  virtual  strain  tensor.  This  principle  obvi¬ 

ously  includes  the  one- dimensional  beam  equations  as  a  special  case;  in  fact  for  a 
virtual  rigid  displacement  of  the  form 

6s  =  6u(s)  —  b  X  64>(s) 


the  virtual  strain  tensor  vanishes  identically  and  from  eqn  IS  we  obtain 

d 

—  {6uT+6<t)-M)  =  0  (19) 

ds 

where 

T  =  y  pdA  M  =  J  b  X  pdA 

are  respectively  the  cross-sectional  stress  resultant  and  moment.  Moreover,  since 
Ss  is  rigid,  we  must  have 


d6u 

ds 


=  64>  xt 


d6d> 

ds 


=  0 


and  substituting  these  expressions  into  eqn  19  we  can  easily  obtain  the  one¬ 
dimensional  beeun  equations 


dM 

ds 


+  txT  = 


0 
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which  show  the  consistency  of  the  three-dimensional  and  one-dimensional  ap¬ 
proaches. 

In  order  to  perform  a  finite  element  approximation  of  eqn  18  it  is  much  easier 
to  rewrite  it  in  the  pulled-back  form, 


j  tr(Si  •  i)  •  e)y/gdA  =  ^  j  Si  ■  pdA  (20) 

in  which  the  hat  denotes  the  pull-back  operation  to  the  cross-section  at  s  =  0,  i.e. 

Si  —  -Ss  p  =  RF  •  p 

Si  =  R^  Se  R  i=^R^  e  R 


and  D  is  the  pulling  back  at  s  =  0  of  the  elastic  tensor,  which  is  constant  with 
resp)ect  to  s  by  our  assumption.  This  is  equivalent  to  refer  all  vectors  and  tensors 
to  the  cross-sectional  components. 

The  finite  element  approximation  of  eqn  20  is  easily  obtained  by  means  of  the 
following  displacement  representation: 

i«‘)  =  jv(r)-s(3) 


where  S  denotes  the  vector  array  of  the  displacement  at  the  nodes  of  the  finite 
element  mesh  of  the  discretized  cross-section,  and  N  denotes  the  matrix  of  the 
shape  functions. 

Using  the  definition  previously  given  of  covariant  derivative  along  the  cross- 
sectional  axes,  the  independent  components  of  the  strain  tensor  e  can  be  given 
the  form  Ai^°)  •  S  {s)  +  B(^“)  •  S{s)  where  JL,  B  are  suitable  given  functions  of 
the  cross-sectional  coordinates  and  ()'  denotes  the  differentiation  along  the  beam 
aixis.  Discretization  of  eqn  20  then  leads  to  the  following  set  of  ordinary  differential 
equations: 


MS'  +  C^S  =  P 
C  s'  +  E-S  =  P 


(21) 


where 

M  =  sa^d-Ma 
C  =  fB^DMA 
E  =  fB^  D  BdA 
P  =  JN'^.pdA. 


The  variables  P  can  be  cadled  nodal  stresses  in  the  sense  that  they  represent  the 
discretized  stresses  acting  upon  the  cross-sectional  surface.  Equation  21  constitutes 
a  set  of  first-order  ordinary  differential  equations  in  the  variables  S  and  P  and 
can  be  put  in  a  second-order  form  by  eliminating  the  nodal  stress  variables.  After 
differentiation  of  the  first  of  eqn  21  and  substitution  into  the  second,  we  obtain 


Equations  21  and  22  axe  the  discretized  differential  homogeneous  equations  of  equi¬ 
librium  with  constant  coefficients  and  symplectic  structure,  hence  the  eigenvalues 
are  of  the  typ>e;  A  =  ^  ib.  From  physical  considerations,  according  to  Saint- 

Venant’s  principle,  we  expect  six  eigenvaJues  with  zero  real  part,  representing  the 
solutions  corresp>onding  to  nonzero  stress  resultants,  and  eigenvalues  with  nonzero 
real  part,  corresponding  to  the  exponentially  decaying  self-equilibrated  modes.  Fol¬ 
lowing  Giavotto  ei  al.  (1983),  the  non-decaying  solutions  are  called  central  solutions 
while  the  decaying  ones  are  csJled  extremity  solutions. 

Considering  the  beam  ais  a  one- dimensional  continuum,  we  are  interested  in  the 
cross-sectional  constitutive  equations,  i.e.  in  finding  out  the  cross-sectional  rigidity 
matrix.  In  this  case,  it  is  common  to  understand  the  strain  energy  as  a  function  of 
the  stress  resultant  and  moment  only,  so  that  we  are  mainly  interested  in  central 
solutions  which  can  be  viewed  as  a  particular  solution  under  a  prescribed  combina¬ 
tion  of  the  stress  resultant  and  moment  at  a  given  cross-section,  zmd  subjected  to 
particular  boundary  conditions  that  avoid  the  existence  of  any  extremity  solution. 
This  goal  can  be  achieved  with  the  procedure  shown  in  the  subsequent  paragraphs. 

Before  continuing,  it  is  extremely  important  to  establish  some  fundamentail  prop¬ 
erties  of  the  operators  M ,  C  and  E.  These  properties  come  from  the  consideration 
that  for  rigid  displacement  of  the  entire  beam,  the  stress  vector  p  is  identically  zero 
so  that  the  nodal  stress  vector  P  is  ailso  zero.  For  our  purposes  the  displacement  i 
can  be  expressed  as  an  infinitesimal  rigid  displacement  of  the  form 

s  =  u(s)  —  b  X  ^{s)  6  =  ■  b 

so  that  the  nodal  displacements  can  be  expressed  as 

S{s)  =  Z-X(s)  X(s)  =  (u(s),<^(s)) 
where  Z  collects  the  constant  operators  Z,  of  the  nodes  Q,. 

Z,  =  [/,-6,  x/]  b.  =  FT  b.. 

Recalling  that  for  rigid  displacement  we  have 

u  =  — c  X  it  —  t  X  (j>  =  ~c  X  0 

or 

-  /  -  T 

X  =  r  X, 

from  eqn  21  it  follows  that  the  operators  M,  C  and  E  must  satisfy  the  properties 

Af  •  z  •  -*-  c’'  •  Z  =  0 

CZT^-|-JFZ  =  0 

and  since  the  operator  M  is  nonsingular  we  also  obtain; 

M~'  ■z  =  -z  r'^ 

.  _i  X  T  -  '“O; 

E  Z=C  M  C  Z 
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2.6  Displacement  Resolution 

Since  it  is  well  known  that  the  central  solutions  entzul  cross-sectioned  distortion,  it 
is  extremely  significant  to  resolve  the  displacement  a  into  a  rigid  part  that  does 
not  strain  the  cross-section  and  into  a  residual  part  w  called  warping,  which  is  re¬ 
sponsible  for  the  cross-sectional  distortion.  This  resolution,  in  pulled-back  notation, 
reads 

s(^*)  =  tt(s)  —  6  X  0(s)  +  w  =  ■  w.  (24) 

It  is  immediately  seen  that  this  resolution  is  six  times  indeterminate,  meaning  that 
the  warping  itself  needs  further  specification.  Obviously,  the  specification  of  the 
warping  entails  the  specification  also  of  the  translation  u  and  the  rotation  be¬ 
cause,  according  to  eqn  24,  they  must  represent  the  same  displacement  s.  Moreover, 
irrespectively  of  the  particular  warping  determination  adopted,  the  stress  distribu¬ 
tion  over  the  cross-section  is  not  affected. 

If  the  same  kind  of  resolution  is  adopted  for  the  virtual  displacement,  the  virtual 
work  principle,  eqn  20  is  rewritten  as: 

J  tr{6e  •  i)  •  €)y/gdA  —  J  ‘ 

+  {Su  +  cx  8u  +  i  X  8^)  ■  T  +  {8^'  -I-  c  x  8^)  •  M.  (25) 

It  can  be  seen  that  in  general,  the  warping  contributes  directly  to  the  internal 
virtual  work,  but,  if  the  interest  is  focused  on  central  solutions  only,  there  is  a  great 
interest  in  defining  the  internal  energj*,  the  warping  u;  and  the  stress  vector  p,  as 
a  function  of  the  stress  resultant  smd  moment  only.  To  this  end,  eqn  25  suggests 
an  intrinsic  definition  of  the  warping  as  a  particular  displacement  ib  satisfying  the 
following  orthogonality  condition: 


j  w  ■  pdA  =  0  (26) 

ais  proposed  by  Borri  and  Merlini  (1986).  By  doing  so,  the  eissociated  quantities 
u  -\-cxu  +  tx<t>  and  0  +  c  x  0,  which  can  be  grouped  in  the  following  vector, 

5  =  (u'  -I-  c  X  u  -t-  f  X  0, 0  -H  c  X  0)  =  X'  -  •  X  (27) 

2issume  the  meaning  of  intrinsic  meaisures  of  the  cross-sectional  deformation,  being 
conjugate  of  the  generalized  stress  resultant  S:  that  becomes  evident  from  eqn  25, 
here  written  in  the  form 

j  tr{8i  ■  D  •  €)y/gdA  =  8£  •  S. 

The  meaning  of  the  orthogonality  condition  in  eqn  26  is  clear:  it  aJlows  one  to  obtain 
the  one-dimensional  expression  of  the  straiin  energy  in  a  straightforwzu-d  manner. 
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Indeed,  it  could  be  observed  that  that  condition  is  too  strong,  since  an  arbitrary 
constant  value  for  that  integrjil  would  suffice:  for  our  purposes,  let  choose  a  null 
value  for  such  constant. 

In  practice,  the  a  priori  imposition  of  the  orthogonality  condition  in  eqn  26  is 
a  difficult  task.  Hence  it  is  better  to  compute  a  particuleir  warping,  which  does 
not  satisfy  eqn  26,  and  then  recover  the  intrinsic  warping  by  means  of  a  special 
projection  as  shown  later  on. 

2.7  Central  Solution 

Let  us  perform  a  finite  element  approximation  of  the  warping  by  writing 

w(«‘) = m°)  • 

where  W  denotes  the  vector  of  the  warping  displacement  at  the  nodes  of  the  finite 
element  mesh.  The  vector  of  the  nodal  displacements  assumes  the  form: 

S  =  Z  X  +  w. 

The  discretized  counterpart  of  the  virtual  work  principle  eqn  25  leads  to  the  follow¬ 
ing  set  of  ordinary  differential  equations: 

K  Y  =  L  (28) 

where 

y  =  {w\w,x',X) 

.  . ,  . r  .  -T  ‘  I 

L  =  {P,P,Z  P,Z  P) 
and  the  matrix  K  hsis  the  following  structure: 

‘  M  MZ  C^Z 

C  E  CZ  EZ 

zi^  M  Z^C^  Z^  MZ  Z^  z 
.  Z^ C  Z^ E  Z^CZ  Z^ EZ  . 

Taking  into  account  the  properties  of  the  operators  M ,  C  smd  E  previously  ob¬ 
tained,  eqn  23,  and  the  definition  of  the  cross-sectional  deformation,  eqn  27,  eqn  28 
transforms  to 

k  Y  =  L 

where 

Y  =  iw\w^e,x) 

£  =  (P,P',5,0) 
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and  the  matrix  K  has  the  following  structure: 

M  M  Z  O 

C  E  C  Z  O 
Z^-M  Z^-C^  Z^  M  Z  O 

o  o  oo 

which  clearly  shows  that  the  rigid  displacement  parameter  X  will  remain  indeter¬ 
minate. 

From  this  equation  the  cross-sectional  deformation  vector  S  is  easily  eliminated, 
leading  to 

M  W'  +  W  =  P  -  M  Z  F~'  S 
C  W'  +  E  W  =  P'  -C  Z  F~^  S 
£  =  F''  ■  (s  -  Z^  ■  M  W'  -  Z^  ■  -W^ 

where: 

F  =  Z  -M-Z 

M  =  M  -  M  ■  Z  ■  F'^  ■  Z^  ■  M 
C  =  C -C  Z  F'^  ■  Z^  M 
E  =  E-C-Z-F~'  -6^. 

Eliminating  P  and  P  by  differentiation  of  the  first  group  and  by  substitution  into 
the  second  group,  tahing  eqn  12  into  account,  the  following  second-order  differential 
equation  set  is  obtained: 

M  -W"  +  {C^  -C)-W'  -E-W  =  G-S  (29) 

where  G  =  C  •  Z  ■  F  -I-  M  ■  Z  •  F  ■'T . 

This  set  of  equations  is  six  times  redundant  since  we  know  that  the  warping  is  six 
times  indeterminate  and  in  fact  it  can  be  recognized  that  the  following  properties 
hold: 

Z  •  Af  =  0  Z  ■  C  =  0 

z^  •  c  =  0  z^  >  f;  =  0 

Z  •  G  =  0. 

Then  it  is  possible  to  eliminate  six  rows  from  the  previous  set  of  equations,  set 
six  components  of  W  ,  W  ,  W  to  zero,  and  solve  the  remaining  equations.  This 
is  equivalent  to  putting  six  independent  single-point  constraints  on  the  warping: 
a  common  practice  to  do  this  in  finite  element  procedures  is  to  add  extremely 
high  zurtificisJ  rigidities  to  the  diagonal  terras  of  six  indepiendent  equations  and  set 
the  corresponding  right-hand  side  to  zero.  The  wau'ping  W  so  determined  will  be 
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denoted  with  an  overbar,  and  the  same  notation  will  be  adopted  for  the  operators 

involved.  ^  _ 

Looking  at  the  form  of  S  given  in  eqn  16,  a  particular  solution  for  W  can  be 
computed  along  with  the  corresponding  cross-sectional  strain  £  and  the  nodal  stress 
vector  P.  This  solution  can  be  concisely  written  in  the  following  form; 

iv  =  w  5  e  =  c  s  p  =  n  s 

where  the  o^rators  W,  C  and  77  are  constant.  Moreover,  since  by  definition  we 
have  S  =  Z  •  P,  the  following  property  holds: 

n  =  I.  (30) 

It  must  be  noted  that  W,  C  strictly  depend  on  the  particular  set  of  single-point 
constraints  imposed  on  the  warping.  On  the  contrary,  the  stress  distribution  is  not 
affected;  consequently,  JT  does  not  depend  on  this  particular  choice.  As  it  will  be 
shown  in  the  next  paragraph,  from  the  knowledge  of  these  qujintities  it  is  p>ossible  to 
determine  an  intrinsic  warping  and  an  intrinsic  cross-sectional  strain,  as  proposed 
by  Borri  and  Merlini  (1986). 

The  particular  solution  previously  mentioned  is  obtained  recalling  the  form  17 
of  the  stress  resultant  and  moment,  and  expression  15  of  the  integration  constants. 
This  form  suggests  one  should  single  out  the  particular  solution  of  the  warping 
displacement  in  the  following  form: 

W  =  Uo  +  <pVo  +  {Uc  -h  <f>Vt)cos<p  +  {U i  -f-  <pV 

Then  from  eqn  29  the  following  sets  of  linear  algebric  equations  in  terms  of 
Vo,  XJc,  Vc,  I/,,  V,  aure  obtained: 

£  •  Vo  =  -G  •  Bo 

EUo  =  -G-Ao  +  cHVo 

{c^M  -I-  B)  -  V,  -H  cH  •  V,  =  -G  •  Be 

_  _  _  _  _  _  (31) 

cH^  ■  Ve  +  (c^M  +  E)-V,  =  -G-Bs 

(c^M  -f- B)  Uc  +  cH  Us  =  -G  A,- cHV,  +  2c^M  ■  V, 

cH^  •  U,  +  (c^M  +  E)-U,  =  -G-  A,  +  2c^M  ■  V,  -|-  cH^  ■  V, 

where  for  sake  of  coinciseness  we  have  defined  H  =  =  C  —  C^.  Recalling 

expressions  15,  these  equations  leawl  to  the  solution  in  terms  of  the  stress  resultant 
amd  moment  corresponding  to  the  section  at  s  =  0,  i.e.  B(0),  and,  in  particular. 
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we  obtain  W(0)  =  Uq  +  and  W  (0)  =  c(Vq  +  +  Us),  that  C2ui  be  used  in 

order  to  compute  the  global  section  strain  f  (0).  The  particular  case  of  stradght 

beams,  i.e.  c  =  0,  can  be  recovered  considering  that  in  this  case  —Vq  =  sVq  and 

Uc  =  Us  =  V,  =  Vs  =  0. 

It  is  extremely  interesting  to  observe  the  symmetry  of  eqns  31  and  that  they  can 
be  solved  in  sequence,  i.e.  at  first  for  Vo  then  for  Uo,  subsequently  for  Vc  and  V, 
and  at  last  for  Uc  and  C/,.  In  addition,  due  to  the  symmetry  of  these  equations, 
we  observe  the  possibility  of  stacking  the  equations  for  Ve,Vs  and  for  Uc,Us  in 
complex  form,  which  can  be  more  efficient  for  computational  purposes.  To  this  end 
let  us  define  the  following  complex  constants: 

U  =  Uc  +  iUs  V  =  V,  +  iVs. 

Then  the  complex  equations  for  U  aind  for  V  are  written  in  the  following  form  that 
replaces  the  last  four  equations  31: 


(c^M  +  icH +  E)-V  = -G- B 

(cW  +  icH  +  E)-U  =  -G-  A  +  c(H  -  2icM)  ■  V. 

2.8  Intrinsic  Warping 

As  intrinsic  warping  w,  according  to  eqn  26,  we  defined  that  particular  warping 
which  is  orthogonal  to  the  stress  distribution.  Then  we  have  P  ■  TV  =  0  for  any 
independent  set  of  stress  resultants  and  moments,  which  entails  il^  •  VV  =  0.  Using 
this  property  allows  us  to  recover  from  the  warping  w  the  intrinsic  wauping  w.  In 
fact,  the  nodal  displacement  vector  S  can  be  written  both  in  terms  of  TV  and  in 
terms  of  the  intrinsic  warping  TV  in  the  following  way: 

5  =  ZX  +  TV  =  ZX  +  TV 


from  which  we  easily  see  that  the  two  warpings  differ  by  a  rigid  displacement  of  the 
cross-section.  In  fact  we  have 

TV-TV  =  -Z.(X-X) 

and  since  both  wairpings  refer  to  the  same  loading  conditions  we  also  have 

TV-TV  =  (W-VV)-5  =  Z-(X-X)  =  Z-  C-  5 

where  Q  is  the  constant  operator  that  takes  into  account  the  dependency  relation 
between  the  warping  distributions; 

W  -  W  =  Z  •  Q. 
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Left-multiplying  this  equation  by  iT^,  we  obtain  the  final  result 

Q  =  /T^  W 

where  the  proi>erty  30  has  been  taken  into  accoimt.  Finally,  the  intrinsic  warping 
distribution  is  expressed  in  the  form 

VV  =  W  -  Z  •  Q  =  “P  •  W 

where  P  =  /  —  Z  •  is  the  particular  projection  needed  to  recover  the  intrinsic 
warping.  That  the  operator  P  is  a  true  projector,  in  fact,  is  easily  verified  since 
P^  =  P.  Moreover,  this  projector  is  orthogonal  to  the  rigid  displacement  of  the 
cross-section  since  we  have  P  •  Z  =  0,  and  its  transpose  is  orthogonal  to  the  stress 
distribution  since  P^  •  i7  =  0. 

After  the  intrinsic  warping  is  determined,  the  intrinsic  cross-sectional  straiin  can 
be  recovered.  This  is  obtained  observing  that 

£  -£  =  X'  -T  -  ■  (X  -X)  =  Q  -  s' -t''  -QS. 

Hence,  tadcing  the  equilibrium  equation  12  into  account  and  setting  £  =  C  •  S,  vfe 
obtain  the  cross-sectional  compliance  matrix  in  the  form: 

C  =  C-Q.t-t^  Q. 

An  alternative  approach  is  also  possible,  eis  indicated  by  Giavotto  et  al.  (1983): 
this  consists  of  computing  the  internal  energy  stored  in  the  section  for  each  inde¬ 
pendent  component  of  the  generalized  stress  resultant  S. 

The  compliance  matrix  C  is  not  singuleu*,  thus  the  cross-section  rigidity  K  is 
obtained  by  inversion. 
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3  Formulation  of  Finite  Elements 


The  Finite  Elements  developed  for  the  cross-sectional  analysis  are  depicted  in  Fig.  3; 
they  are  isoparametric  elements  on  the  cross-section  domain.  The  most  general 
element  is  the  Plane  Element:  it  is  a  full  three-dimensional  solid  element,  whose 
view  on  the  cross-section  is  a  quadrilateral  polygon;  the  Lamina  Element  is  quite 
similar  to  it  and  is  designed  to  model  very  thin  solids.  The  Panel  Element  is  a  plane 
stress  element  whose  view  on  the  cross-section  is  an  arc  of  line:  it  is  used  to  model 
walls  in  hollow  beams;  the  Joint  Element  is  quite  similar  to  it  and  only  works  in 
shear.  Finzdly,  the  Stringer  Element  is  a  one-dimensional  element  which  is  seen  on 
the  cross-section  as  a  point. 

3.1  Plane  Element 

This  element  is  a  quadrilateral  element  used  to  model  the  cross-section  of  solid 
beams.  The  formulation  of  this  element  is  based  on  an  analytical  expression  of  the 
strain  and  stress  tensors  in  terms  of  cross-sectional  components,  that  is  given  in  the 
following. 

From  the  programmer's  point  of  view  it  is  easier  to  adopt  a  little  different  nota¬ 
tion  than  that  used  so  far:  in  particulau:  let  us  denote  with  {V}  the  vector  atrray  of 
the  cross  sectional  components  of  a  generic  vector  V,  i.e.  {V}^  =  [Vi,  ^3]  where 
14  =  e*  •  V,  Following  the  same  kind  of  notation,  we  denote  with  Cik  =  £{  •£•  the 
cross-sectional  components  of  the  strain  tensor  c,  and  for  convenience  we  introduce 
the  strain  array  vector  where 

{^o}^  =  [2^13'  2e23i  C33]  =  [cii?  €221 -^u] 

represent  the  out-of-plane  and  in-plane  strain  tensor  arrays. 

Taking  into  account  the  expressions  of  the  displacement  gradient  tensor  reported 
above,  we  can  write 

and 

{€,}  =  I'D, US) 

where  {5}  is  the  displacement  vector  array  and  the  differential  operators  [T>o\  and 
[Di]  have  the  following  definition: 

"P  -C3  cj  + 

^3  V  -C,  -h 

-C2  Cl  T> 
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ELEMENT  LIBRARY 


NAME 

TOPOLOGY 

PLANE  ELEMENT 

2 

3  Aw  2  H  — ^ 

/ 

10^^ 

1 

LAMINA  ELEMENT 

i 

1 

3 

STRINGER  ELEMENT  ! 

j 

•  1 

JOINT  ELEMENT 

1  •-•  2 

THIN  PANEL  ELEMENT 

! 

1  •  ■  •  2 

Figure  3:  Element  library. 
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with  "D  defined  as  follows: 


The  metric  determinant  ^/g  can  be  computed  by 

y/g  =  h-  CjC  + 

and  we  recall  that,  within  the  approximation  of  a  helicoidal  be2an,  the  vector  arrays 
{i}  «ind  {c}  are  constants. 

Assuming  the  following  finite  element  discretization 

{s(e)}  =  mm  ■  isis)} 


and  letting  S  =  -^-r,  we  can  write 

£74'* 


{^o}  =  }  4-  [So]{5} 


{£,)  =  18, HS), 


where: 

i>toi  =  ;4iMr)i  isoi  =  ■^[i>oiiJV(r)i 

Collecting  the  previous  formulae  we  obtzdn 


(8,1  =  |i>,][jv(ni. 


{i}  =  (.4|{S  )  +  18]{S} 


where: 


[^1  = 


[B]  = 


The  numerical  evaluation  of  the  matrices  [Bq]  and  [B/]  requires  the  differentiation 
of  the  shape  funtions  [iV(4“)]  with  respect  to  the  cross-sectional  coordinates  (f“: 
this  is  easily  performed  according  to  the  particular  type  of  finite  element  used. 
The  element  implementation  follows  an  isoparametric  formulation  beised  upon  2 
to  4  nodes  for  each  side  and  aJlows  for  linear,  parabolic  and  cubic  displacement 
interpolation.  The  toted  number  of  nodes  remges  from  4  to  12.  The  shape  functions 
eure  rep>orted  in  Table  1,  where  p  emd  q  denote  non-dimensional  cartesian  coordinates 
on  the  parent  domeun,  as  shown  in  Fig.  4. 

The  constitutive  equation  a  =  D  •  e  cem  be  written  in  terms  of  out-of-plane  and 
in-plane  eomp>onents,  with  reference  to  the  cross-section  triad,  in  the  following  way: 

{^o}  =  l^oo]{co}  4-  [£>o/]{c/} 

{&i)  —  [Dio]{eo)  4-  [!)//]{£/} 
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Element  order 


Shape  function 


JV,  =  }(l+p)(l+8)  +  KJVs  +  iV,) 

=  JO-p)(i+«)  +  ;W  +  ^e) 

%=  J(l-p)(l-«)  +  |(iV,  +  JV,) 
-*<  =  J(i+p)(i-«)  +  |(Wt  +  %) 
Ns  =  l(l-p')(l  +  ,) 

^6  =  i(l  -p)(l 
Arr  =  i(l-p“)(l-9) 

JV,=  l(l+p)(l-9=) 


JV,  =  i(l  +p)(l  +  9)(-10  +  9(p'  +  q^)) 
N,  =  i(l  -  p)(l  +  5)(-10  +  9(p"  +  9*)) 
'V3  =  ^(1  -  P){1  -  9)(-10  +  9(p2  +  9*)) 
=  »('  +PK>  -9)(-10  +  9(p^  +9*)) 
iVs  =  ^(l+9)(l-p')(l  +  3p) 

JV8  =  A(l  +  9)(l-p")(l_3p) 

.IV,  =  ^(l-p)(l-9')(l  +  3p) 
JV.=  ^(I-pK1-9")(1-3p) 
^9=i(l-9)(l-P*)a-3p) 

Pf.o  =  5(l-9)(l-p")(l+3p) 


W„  =  ^(l+p)(l-9')(l-3p) 
^>2  =  P)(^  “  9^)(1  +  3p) 


Table  1:  Shape  functions  of  the  Plane  Element. 
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Figure  4;  Parent  domain  of  the  Plame  Element. 

The  array  vectors  {&o}  *  ^md  {o’!}  are  the  out-of-plane  and  in-plane  stress  ten¬ 
sor  arrays,  while  the  matrices  [.Doo]i  [Doi]-,  [D[o],[Dii],  group  the  cross-sectional 
components  of  the  elasticity  tensor  D.  Since  D  is  symmetric,  the  equivalence 
[Doi\  =  [DioY  is  assumed.  These  matrices  depend  on  the  particular  material: 
the  user  may  supply  the  elasticity  tensor  D  directly  or  its  inverse,  or  he  can  sup¬ 
ply  the  elastic  properties  along  the  orthotropy  directions  for  orthotropic  materials. 
Orthotropic  materials  are  allowed  with  an  orthotropy  plane  normal  to  the  cross- 
section:  the  user  must  supply  the  orientation  of  this  plane  on  the  cross-section  and 
the  orientation  of  the  orthotropy  directions  in  this  plane  by  means  of  two  angles. 
These  angles  are  unique  for  each  element  and  can  be  supplied  with  reference  to 
either  the  physical  or  the  p2U'ent  domain. 

The  standard  integration  procedure  uses  Gauss-Legendre  points  and  weights  to 
compute  the  stiffness  matrices  [Af],  [CJ  arid  [E]  of  the  element.  In  the  case  of  linear 
and  peirabolic  interpolation  functions  it  is  possible  to  force  the  use  of  integration 
points  located  at  the  nodes  of  the  element:  this  allows  us  recovering  stresses  at  the 
nodes  themselves. 


3.2  Lamina  Element 

This  element  is  very  similar  to  the  plane  element  and  has  been  designed  to  model 
very  thin  structural  elements,  as  for  example  interlaminar  layers.  It  can  be  consid¬ 
ered  as  a  degeneration  of  the  plane  element,  in  fact  the  formulation  is  almost  the 
same,  except  the  strain  which  is  kept  constant  through  the  element  thickness.  The 
shape  functions  along  the  lamina  direction  are  the  same  as  the  plane  element,  while 
a  linear  interpolation  is  used  in  the  transverse  direction.  A  relaxed  integration  can 
be  used  through  the  thickness.  The  geometry  of  the  lamina  element  is  defined  by 

^The  array  vector  {&o  }  representa  the  stress  vector  acting  at  the  generic  point  on  the  cross-section 
and  coincides  with  the  vector  elsewhere  called  p. 
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mezms  of  2  to  4  nodes  on  its  reference  line  and  by  means  of  its  thickness  at  each 
node.  This  element  can  only  be  orthotropic  with  an  orthotropy  plane  tangent  to 
its  reference  line. 


3.3  Panel  Element 

The  panel  element  can  be  used  to  model  thin- walled  beams:  it  is  a  thin  element 
of  constant  thickness,  defined  from  a  geometrical  point  of  view  zis  a  surface  which 
intersects  the  cross  section  along  a  line  /(^“).  An  isoparametric  description,  with 
a  number  of  nodes  ranging  from  2  to  4,  is  assumed  for  this  line. 

From  a  mechanical  point  of  view  the  panel  is  a  membrane  working  in  plane 
stress;  then  letting  v  be  the  unit  normal  to  the  line  the  identity  •  i/  =  0 

is  trivizdly  satisfied.  Moreover,  letting  be  the  unit  tangent  to  the  line  /(^"),  the 
significant  components  of  the  stress  tensor  can  be  grouped  in  the  stress  zuray  vector 

{arp]^  =  ^m3>33] 

being  =  fi  a  ■  Cs  zind  ^33  =  63  •  6^  •  63.  Following  the  same  kind 

of  notation,  we  define  the  strain  tensor  array  {cp}  zis 

{ep}^  =  £33] 

and  the  plane  stress  state  constitutive  equation  as 

{&p}  -  [bp]  ■  {6p} 

where  [Dp]  denotes  the  elzisticity  matrix  typical  of  the  plane  stress  state. 

The  expression  of  the  components  26,^3  zmd  C33  comes  from  the  general 
expression  of  the  strain  tensor  taking  into  account  that  in  this  case  we  have 


where  ds^  represents  the  differentizd  of  the  curvilinear  abscissa  along  the  line  /. 
Then  the  paxtizil  derivatives  with  respect  to  the  cross  section  coordinates  are  to  be 
replaced  by 

a  a 
ae 

and  from  the  general  expression  of  the  displacement  gradient  tensor  we  obtain: 


I  (  ds  ds\ 


Assviming  the  following  finite  element  discretization 

(»)  =  (iV(^.)l  ■  («(»)) 
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Element  order 

Shape  function 

Linear 

=  ^(i-p) 

Quadratic 

N,  =  ‘(p*  -  p) 

N,  =  i(p'  +  P) 

A  =  (l-P=) 

Cubic 

jV,  =  ±(-l-hp(l+p(9-9p))) 

■^2  =  +  P(9  +  9p))) 

iVa  =  ^(9  +p(-27  +  p(-9  -h  27p))) 
N,  =  i(9  -h p(+27  +  p{-9-  27p))) 

Table  2:  Shape  functions  of  the  Pzmel  Element. 


we  can  write 

{£>}  =  [Xp1[jv(.„)Hs'}  +  |B?)i;v(s.)j{s} 
where,  as  it  is  easily  verified, 


0  0 

h-  h. 

y/9  y/9 

0  0 


0 

0 

1 


^r9\ 


[Bp]  = 


ds. 


M2 


ds.. 


M2C3  -  ani 


-C2 


ds^ 


-M1C3  -  an2 


Cl 


ds. 


M1C2  ~  M2C1  + 


ds^ 


—a 


ds. 


fx  • 

and  a  =  — fia  =  ■  Cq.  The  shape  functions  for  linear,  parabolic  and  cubic 

.  ^ 

interpolation  are  reported  in  Table  2,  where  p  denotes  non-dimensional  cartesian 
coordinates  on  the  parent  domain,  as  shown  in  Fig.  5. 
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X, 


Figure  5:  Parent  domain  of  the  Panel  Element. 

When  using  a  curvilinear  element,  care  must  be  taken  for  two  reasons.  First,  the 
matrix  [Dp]  has  constant  components  referred  to  the  local  directions  /i,  v  and  63, 
hence  it  is  not  a  constant  matrix  when  referred  to  the  cross  section  triad.  Second, 
due  to  the  component  this  element  requires  different  from  zero  in  order 
to  be  equilibrated:  this  stress  component  can  be  provided  by  other  elements.  The 
program  takes  care  of  the  first  problem:  the  user  is  allowed  to  supply  the  elasticity 
matrix  [jD/.]  with  reference  to  the  parent  domain.  About  the  second  problem,  it  is 
left  to  the  user  to  judge  the  opportunity  of  using  curvilinear  elements.  Nevertheless, 
it  is  always  possible  to  supply  singular  data  for  [Dp]  so  that  will  be  zero. 
Furthermore,  in  order  to  emulate  semi-monocoque  behavior  the  user  can  set  to  zero 
both  the  extensional  rigidities. 

The  panel  element  can  be  also  used  in  order  to  model  a  thin  laminate:  its 
stiffness  matrix  can  be  evaluated  stacking  the  characteristics  of  a  sequence  of  layers 
that  can  be  made  each  of  different  material  with  different  orthotropy  orientations. 
Obviously  this  element  can  be  used  only  if  the  interlaminar  stresses  are  considered 
not  to  be  important. 

3.4  Joint  Element 

The  joint  element  has  been  implemented  to  model  an  elastic  link  between  an  offset 
stringer  and  a  thin  panel.  It  is  a  two-node  shear  element  obtained  as  a  degeneration 
of  a  panel  element,  with  both  the  extensional  rigidities  set  to  zero. 

3.5  Stringer  Element 

The  stringer  element  is  mainly  used  to  model  stiffners  in  thin-walled  beams;  on 
the  cross-section  domain,  this  element  can  be  classified  as  a  point  element.  A 
stringer  is  modeled  by  a  line  element  which  is  able  to  sustain  an  axial  load  N  along 
the  direction  A  of  the  stringer  itself  with  an  axial  rigidity  K  .  Letting  sa  be  the 
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curvilinear  abscissa  along  the  stringer,  the  constitutive  equation  of  the  stringer  is 
obviously  given  by 

iV  =  A'A0A~ 

dsx 

where  a  is  the  diplacement  of  the  stringer. 

In  order  to  use  this  element  in  the  context  of  the  present  approach,  we  have  to 
specify  the  derivative  of  the  displacement  with  respect  to  the  curvilinear  abscissa 
in  terms  of  the  derivative  along  the  beam  axis  and  evaluate  the  virtual  change 
of  the  potential  energy  of  the  stringer.  To  this  end,  let  Q  be  the  point  of  the 
intersection  of  the  stringer  with  the  generic  cross-section  and  let  be  the  cross- 
sectioned  coordinates  of  the  point  Q.  Moreover,  if  dQ  denotes  an  infinitesimal  length 
vector  along  the  stringer,  we  have 


dQ  =  g^di’‘  =  Ads  A, 


(standing  on  the  cross-section,  this  element  can  be  classified  as  point  eleTnent)and 
obviously  dsx  =  \/dQ  ■  dQ.  Along  the  line  element  of  the  stringer  the  coordinate 
differentials  d^*  can  be  expressed  in  terms  of  dsx  as  follows, 

d^*  =  (7^  •  Xdsx, 

and  assuming  the  stringer  intersect  the  cross-section  (i.e.  A  •  C3  different  from  zero) 
we  obtain  the  relationship  between  d(^  eind  dsx  in  the  form 

d^^  =  ■  Xdsx  =  J^ds.x, 

n/^ 


whence 


da  _  e^  -  X  da 
dsx  y/g  d^3- 


Then,  the  virtual  cheinge  of  the  potential  energy  of  an  infinitesimal  length  of 
stringer  is  given  by: 


^  .  Ndsx  =  ^  •  TVdf  =  A'^A  ®  A  •  ^ . 
dsx  d^3  ^ 


This  equation  shows  that  the  contribution  of  the  stringer  to  the  virtual  change 
of  the  potential  energy  depends  on  the  position  of  the  stringer  itself  within  the 
cross-section,  due  to  the  different  lengths  of  d!^  and  dsx. 

In  terms  of  cross-sectional  components,  since  the  shape  function  reduces  in  this 
case  to  the  identity,  we  have  «  =  5  and  therefore 


<V  =  ({5')  +  [cx  /|{S))’'|M|  ((S’)  +  [c  X  /){S})  , 
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where  [Af]  is  given  by: 


(m1  =  a-5^{a)(a}’-. 


The  other  finite  element  stiffness  matrices  [C]  and  [£?]  are  then  computed  by  the 
formulae: 

[C]  =  (c  X  lf[M]  (£?]  =  [c  X  lf[M][c  X  /]. 
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4  The  Computer  Code 

This  section  draws  the  basic  ideas  of  the  orgamization  of  the  computer  code  as 
it  should  be  for  an  efficient  analysis  of  the  cross-section  of  curved  and  twisted 
anisotropic  beams.  The  numerical  examples  discussed  later  on  have  been  obtained 
by  an  early  program  written  to  analyze  only  straight  beams  then  modified  for  the 
present  scope.  The  modifications  performed,  even  if  correct,  are  not  the  best  from 
the  point  of  view  of  program  architecture  and  modularity:  in  particular,  in  the 
present  version  the  displacement  resolution  is  performed  at  element  level  and  should 
be  brought  outside,  following  the  formulation  given  above.  Thus,  most  but  not  all 
the  features  discussed  herein  are  actually  implemented  in  the  present  code. 

4.1  Input  Data 

Input  data  au'e  grouped  into  three  main  blocks: 

•  Control  data: 

—  caise  description; 

—  execution  flags  and  parameters,  curvature  data; 

—  outputs  required  and  node/element  selection; 

—  eigenvalue  solution  parameters. 

•  Node  data: 

—  node  identifier  and  coordinates: 

—  nodal  renumbering  sequence; 

—  single  point  constraints: 

—  multiple  point  constraints. 

•  Element  data: 

—  material  properties; 

—  element  properties; 

—  element  connections. 

4.2  Disk  and  Memory  Data  Management  System 

The  computer  program  interfaces  with  mass  storage  units  by  means  of  a  home-made 
disk  manager  based  on  a  paging  technique  similar  to  that  used  by  virtual  operating 
systems.  Only  few  physical  files  are  required  by  the  disk  manager:  that  has  the  nice 
effect  to  allow  the  user  handling  simple  procedures  to  perform  start,  restau't,  backup 
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aaid  restore  operations.  On  the  contrary,  several  logical  files  axe  allowed  within  the 
program  and  can  be  accessed  both  in  sequential  and  in  direct  mode. 

A  memory  manager  is  also  implemented  sJlowing  memory  dynamic  allocation, 
memory  packing  and  other  features.  This  memory  manager  is  fully  integrated 
with  the  disk  manager  to  get  a  high  programming  flexibility.  The  performances  of 
this  system  can  be  tuned  at  installation  time,  according  to  the  specific  hardware 
configuration,  by  means  of  a  set  of  parameters  grouped  in  a  block  data. 

4.3  Organization  of  the  Computer  Code 

The  program  is  organized  in  several  main  modules: 

•  Cast  Input.  This  module  reads  the  data  for  the  execution  control: 

-  Case  description,  execution  flags,  parameter  values  and  restart  options; 

—  Outputs  required:  echo  of  input  data,  displacements,  stresses,  strains, 
intermediate  results; 

—  Node/element  data  selection  for  warping,  stress  and  strain  printouts; 

—  Reference  systems  for  input  data. 

•  Data  Input  This  module  performs  the  following  items: 

—  Node  data  input.  The  nodal  data  aure  rearranged  in  ascending  order 
with  regard  to  the  node  identifier;  index  vectors  for  external/internal 
numbering  correspondence  are  prepared,  accounting  for  possible  externail 
reordering  sequence. 

-  Constraint  input.  Single  and  multiple  point  constraints  are  available.  An 
automatic  procedure  is  implemented  for  eliminating  any  singularity  due 
to  rigid  body  motions;  a  table  of  constrained  DOFs  is  printed.  Single 
point  constraints  are  handled  eis  very  strong  springs;  the  default  spring 
rigidity  can  be  modified  in  input.  Each  multiple  point  constraint  generate 
a  kinematic  equation  that  is  added  to  the  analysis  by  means  of  a  Lagrange 
Multiplier. 

—  Equations  numbering.  The  equations  related  to  each  node  are  numbered, 
accounting  for  the  presence  of  multi-point  constraints;  each  constraint 
equation  is  automatically  inserted  eifter  the  first  degree  of  fredoom  in¬ 
volved  in  order  to  ensure  any  singularity  to  be  removed.  The  equations 
related  to  the  generalized  cross-section  deformation  parameters  are  num¬ 
bered  as  last  ones. 

—  Element  Data  input.  The  element  data  are  stored  in  a  logical  file  struc¬ 
tured  with  variable  length  records;  each  record  contains  the  element 
identifier,  the  element  type,  the  connected  nodes,  the  material  data  and 
physical  properties  and  other  informations  specific  to  the  element. 
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—  Equation  pointer  2aTay.  For  each  element  an  array  of  pointers  to  the 
equations  involved  is  set  and  stored  on  disk. 

•  Element  Matrix  Computation.  The  matrices  M,  C,  E  and  P  are  computed 
for  each  element  and  stored  on  a  disk  file. 

•  Assembling  Module.  The  matrices  of  the  cross-section  finite  element  zaialysis 
are  assembled  by  this  module.  The  assembling  process  can  be  performed 
either  with  an  in-core  procedure  or  with  an  out-of-memory  one.  The  out-of¬ 
memory  procedure  runs  with  a  frontal  technique  which  is  very  effective  for 
large  problems.  The  in-core  procedure  is  based  on  a  skyline  assembler/solver 
routine  which  can  speed  up  the  execution  when  large  amount  of  core  memory 
is  available.  The  assembler /solver  package  can  work  either  in  real  or  complex 
mode,  for  the  analysis  of  straight  or  curved  beams  respectively. 

•  Solution  Module.  This  module  performs  the  solution  iterations  according  to 
the  execution  options:  two  solution  steps  on  a  real  matrix  are  run  for  straight 
beams,  while  six  steps  axe  required  on  a  Hermitian  matrix  for  curved  beams. 
In  any  case  the  coefficient  matrix  is  the  same  for  each  step,  thus  only  one 
decomposition  is  needed.  The  right-hand  side  vector  is  setup  at  each  step, 
then  the  frontal  or  skyline  solution  is  performed.  Both  the  implemented  solvers 
contain  an  automatic  procedure  for  singularity  elimination  that  cam  be  turned 
off  with  an  input  option.  Finally,  the  section  intrinsic  warping  is  recovered, 
as  discussed  in  the  theoretical  section  of  this  report. 

•  Cross-section  Stiffness  and  Mass  Module.  The  6x6  stiffness  matrix  of  the 
cross-section  is  computed  and  it  is  evaluated  also  with  reference  to  the  prin¬ 
cipal!  axes  of  the  normal  stress.  By  integrating  the  density  on  the  area  of 
each  element  the  maiss  matrix  of  the  cross-section  is  also  computed.  Then  the 
center  of  gravity  of  the  cross-section  and  its  inertial  properties  and  principal 
axes  axe  computed. 

•  Stress  and  Strain  Recovery.  Stresses  aind  strains  can  be  recovered  either  with 
reference  to  the  element  fraimes  or  to  the  cross-section  global  frame. 

•  Extremity  Solution.  The  solution  of  the  extremity  problem  cam  be  permormed 
with  two  different  formulations:  a  full  amalysis  to  extract  all  the  eigenvaJues 
and  a  reduced  one,  based  on  the  Lanczo’s  algorithm,  to  extract  only  a  few 
eigenvalues  of  interest.  This  module  computes  eigenvalues  amd  eigenvectors 
and  can  recover  the  eigenstresses  too. 

•  Printout  Module.  It  performs  any  output  selected  by  the  user  and  prints  the 
cross-section  global  charaw:teristics. 
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Figure  6:  Elementary  study  cases.  Simply  curved  and  twisted  beams. 

5  Numerical  ExLauples 

5.1  A  Curved  or  Twisted  Isotropic  Solid  Beam 

The  first  study  case  presented  is  very  simple  and  aims  to  only  determine  the  fea¬ 
sibility  of  the  method  presented  without  any  claim  to  be  exhaustive.  An  isotropic 
square  cross-section  of  120x120  mm  is  analyzed  and  results  Eire  presented  for  two 
cases,  see  Fig.  6:  a  simply  curved  beam  with  radius  of  curvature  of  200  mm  (Case 
A)  and  a  simply  twisted  beaim  with  a  twist  of  90°  degrees  over  314.5  mm  of  length 
(Case  B). 

The  beam  cross-section  has  been  modeled  with  25  square  8-node  isoparametric 
elements  for  a  total  of  294  DOFs.  The  analysis  by  the  present  approach  is  com¬ 
pared  with  theoretical  results,  when  available,  and  with  three-dimensional  numerical 
analyses  by  MSC-NASTRAN.  The  three-dimensional  discretization  of  half  a  beam 
is  composed  of  10  or  12  slices  of  25  20-node  isoparametric  solid  elements  (CHEXA), 
for  a  total  of  250  elements  and  4248  DOFs  for  case  A,  and  300  elements  and  5040 
DOFs  for  case  B.  These  models,  see  Fig.  7,  are  loaded  at  the  free  extremity  with 
suitable  nodal  forces  giving  the  desired  resultants.  Note  that  the  present  approach 
involves  more  than  one  order  of  magnitude  reduction  of  the  problem  size  relative 
to  MSC-NASTRAN. 


Figure  7:  Elementary  study  cases.  MSC-NASTRAN  mesh. 

Case  A:  Curved  Beam 

Theoretical  results  cire  only  available  for  a  simply  curved  beam  acted  upon  by  a 
bending  moment  in  the  curvature  plane.  A  comparison  of  the  normal  stress  is  given 
in  Fig.  S:  it  can  be  seen  that  the  present  analysis  and  the  three-dimensional  one 
agree  very  well  and  discrepancies  are  confined  within  0.3%.  We  also  observe  that 
the  maximum  value  of  the  compression  stress  is  about  30%  greater  them  the  straight 
beam  with  the  same  cross-section. 

Figures  9, 10  and  11  report  some  stress  distributions  for  different  load  conditions, 
showing  significzmt  effects  due  to  the  cur\'ature.  .\s  expected,  the  stress  distributions 
are  different  from  the  corresponding  ones  of  the  straight-beam  ceise:  for  example, 
the  in- plane  stress  component  <ru  is  different  from  zero,  the  sheair  stress  distribution 
<ri3  due  to  is  not  symmetric  and  the  shear  center  is  closer  to  the  curvature  center. 
.4  comparison  \uth  MSC-NASTR.4N  results  is  given  in  Fig.  12:  the  difference  in 
computed  stresses  is  under  2%. 

Case  B:  Twisted  Beam 

Also  in  the  twisted  beaun  case  is  the  difference  with  respect  to  the  rectilinear  beam 
with  the  same  cross-section  very  significant.  In  Fig.  13  the  shear  stress  under  four 
elementary  load  conditions  is  plotted:  sheau-  stresses  aurise  at  the  four  corners  even 
under  bending  amd  tension  loading,  thus  satisfying  the  boundary  conditions  on  the 
lateral  surface. 

For  this  case  only  a  comparison  of  the  present  analysis  vs  a  three-dimensional 
one  is  reported,  as  theoretical  results  were  not  found:  the  stress  distributions  on  the 
crpss-section,  see  Fig.  14,  are  quite  similau:  for  the  two  anadyses  even  if  discrepancies 
au-e  more  evident  than  for  caise  A.  A  more  detailed  discretization  would  reduce  these 
differences. 


Figure  8:  Case  A.  Normal  stress  due  to  bending. 


Figure  9:  Case  A.  Stress  <7u  contour  plot  due  to  Mi. 


Figure  10:  Case  A.  Stress  <Ti3  contour  plot  due  to  Tj. 


Figure  11:  Case  A.  Shear  stresses  due  to:  (a)  Ti,  (6)  Ta,  (c)  Mz. 


c  d 

Figure  13:  Case  B.  Shear  stresses  due  to:  (a)  Mi,  (6)  M3,  (c)  Ti,  (d)  Ta. 
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Figure  14;  Case  B.  Some  stress  distributions  due  to  different  section  resultants. 
Comparison  with  MSC-NASTRAN  (dotted  lines)  results. 


5.2  A  Rectilinear  Composite  Box-Beam 

As  a  second  example,  the  rectilinear  composite  box-beam  already  investigated 
by  Smith  juid  Chopra  (1990)  is  considered  (see  Fig.  15):  the  experimented  and 
numerical  results  they  present  give  the  opportunity  to  compare  the  performances  of 
different  approaches  to  the  emalysis  of  comf>osite  beams.  Ten  specimens  have  been 
analyzed  having  two  different  aspect  ratios  and  several  laminations:  the  geometric 
data,  the  stacking  sequences  and  the  elastic  properties  of  the  material,  unidirectional 
GR/EP  AS4/3501-6,  are  reported  in  Tables  3,  4  and  5.  The  stacking  sequences 


Case  A 

Case  B 

L 

30.  in. 

30.  in. 

c 

2.060  in. 

.0953  in. 

d 

1.025  in. 

.0537  in. 

Lfd 

29 

56 

ply  thickness 

.005  in. 

.005  in. 

wall  thickness  (6  plies) 

.030  in. 

.030  in. 

Table  3:  Box-beam.  Geometric  data. 
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Figure  15;  The  rectilinear  composite  box-beam. 


Specimen 

# 

Configuration 

Upper  and 
lower  walls 

Side 

walls 

■BH 

symmetrical 

0/90 

anti-symm. 

15 

15/15 

anti-symm. 

0/30 

anti-symm. 

0/45 

isa 

5 

symmetrical 

15 

15/-15 

6 

symmetrical 

30 

symmetrical 

45 

45/-45 

8 

B 

symmetrical 

15 

15/-15 

9 

B 

symmetrical 

30 

30/-30 

10 

B 

symmetrical 

45 

45/-45 

Table  4:  Box-beam.  Stacking  sequences. 


Orthotropic  Young’s  modulus 

El 

—  20.59  msi 

Transversal  Young’s  modulus 

Ej 

=  1.42  msi 

Tangentizd  modulus 

Glt 

=  0.87  msi 

Poisson’s  ratio 

vlt 

=  0.42 

Table  5:  Box-beam.  Elastic  proi>erties. 
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Figure  16:  Box-beam.  Symmetric  and  Jinti-symmetric  configurations. 


Figure  17:  Box- beam.  Mesh  of  the  cross-section. 
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Figure  18:  Box-beam,  Specimen  3.  Deformed  cross-section  due  to:  (a)  bending 
moment  M],  (b)  torque  M3  and  (c)  axial  force  T3. 


considered  allow  us  to  auialyze  either  uncoupled  behavior  or  symmetric  or  anti¬ 
symmetric  couplings:  refer  to  Fig.  16  for  the  orientation  of  the  laminations.  The 
model  of  the  cross-section,  see  Fig.  17,  stands  on  288  8-node  isoparametric  plane 
elements,  each  bounded  within  a  single  ply:  the  analysis  DOFs  amount  to  2886. 

The  deformed  cross-sections  under  a  few  elementary  load  conditions  are  shown 
in  Fig.  18.  Some  diagonal  terms  of  the  flexibility  matrix  of  the  cross-section,  namely 
flexure  about  Xi  and  torsion,  are  given  in  Table  6  for  adl  the  specimens  and  plotted 
in  Fig.  19  versus  the  characteristic  angle  of  lamination. 

Some  coupling  effects  are  given  in  Table  7:  two  off-diagonal  terms  of  the  flexibil¬ 
ity  matrix  are  rep)orted  for  2dl  the  specimens.  The  section  torsion  under  axial  load 
and  bending  is  also  plotted  in  Fig.  20  versus  the  characteristic  angle  of  lamination. 
These  plots  show  that,  as  it  is  expected,  for  the  symmetric  specimens  (1  and  5 
to  10)  axial  extension  does  not  couple  with  torsion  while  flexure  does;  the  opposite 


46 


1.042  X  10-* 

1.406  X  10'® 

2.665  X  10-* 

1.125  X  10-* 

1.844  X  10-* 

6.445  X  10-® 

2.018  X  10-* 

6.836  X  10-® 

2.890  X  10~® 

9.790  X  10-® 

6.860  X  10-® 

7.260  X  10-® 

1.250  X  10"® 

6.960  X  10-® 

2.380  X  10"® 

8.332  X  10-® 

5.620  X  10*® 

6.120  X  10“® 

1.060  X  lO-** 

5.903  X  10-® 

Table  6:  Box-beam.  Flexure  and  torsion  due  to  unit  moments. 


Table  7:  Box-beam.  Some  coupling  effects. 


Figure  19:  Box-beam.  Flexure  and  torsion  due  to  unit  moments. 
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Figure  20:  Box-beam.  Some  couplings  effects. 
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hapi)ens  for  the  anti-symmetric  specimens  (2  to  4).  Moreover,  the  influence  on  these 
couplings  of  the  fiber  orientation  and  of  the  specimen  2ispect  ratio  can  be  observed. 
It  is  well  known  that  the  knowledge  of  results  like  these  ones  would  be  a  valuable 
support  while  tailoring  a  real  beam,  for  example  for  aeroelastic  purjwses;  all  the 
elzistic  effects  must  be  taken  into  account  in  order  to  achieve  a  desired  coupling 
level,  and  this  can  be  a  very  difficiilt  task  for  a  complex  beam  if  adequate  tools  are 
not  at  hand,  as  this  quite  simple  example  of  a  box-beam  can  demonstrate. 

The  cross-section  stiffness  has  been  integrated  along  the  beam  span  to  obtain 
a  beam  model  for  the  experiments  reported  by  Smith  and  Chopra  (1990);  that 
allows  us  to  extend  the  comparison  among  different  analysis  methods,  presented  in 
that  work.  Figure  21  shows  the  results  obtained;  the  effectiveness  of  the  present 
approach  with  respect  to  others  methods  can  be  appreciated  even  if  we  cannot  point 
out  a  best  one. 

Hereafter  the  strain  and  stress  distributions  across  the  wall  thickness  are  pre¬ 
sented  for  two  points  of  the  cross-section  of  specimen  3,  chosen  at  the  mid-side  of  the 
upper  and  right  walls,  for  two  load  conditions,  namely  axial  load  and  torque.  As  it 
will  be  observed,  some  stress  or  strain  components  are  discontinuous  between  lami¬ 
nae;  accounting  for  these  discontinuities  is  implicit  in  the  present  formulation  based 
on  a  fully  three-dimensional  stress  state.  Obviously  discontinuous  elastic  proper¬ 
ties  must  be  reproduced  with  an  adequate  finite  element  model  of  the  composite 
cross-section. 

In  the  case  of  axial  load.  Figs.  22  to  25,  the  stress  component  cr33  is  discontinuous 
between  layers  at  0®  and  layers  at  30®,  while  strain  €33  keeps  constant.  As  a  conse¬ 
quence,  alternate  shear  stresses  <ri3  in  the  upper-side  wall  and  a’23  in  the  right-side 
one  arise;  the  corresponding  strains  613  2md  C23  are  continuous  and  linearly  varying 
across  the  thickness.  A  transverse  stress  state  can  be  observed;  the  extensional 
stress  <rii  in  the  upper-side  wedl  and  <722  in  the  right-side  one  are  alternate,  while 
the  corresponding  strzuns  en  and  €22  are  linearly  varying.  Transverse  strains  £22  or 
£11  are  present,  while  stresses  <722  or  crn  are  zero.  The  transverse  shear  stress  eri2 
and  strain  £12  are  zero. 

In  the  case  of  torque.  Figs.  26  to  29,  the  results  are  qualitatively  similar;  the 
stress  components  oria  in  the  upper-side  wall  and  <r23  in  the  right-side  one  are 
discontinuous  between  the  layers,  while  the  corresponding  strains  £13  and  £23  are 
linearly  varying,  and  an  axial  stress  <733  arises,  with  null  strain  £33.  A  discontinuous 
behavior  of  the  extensional  stress  Vw  in  the  upi)er-side  wadi  and  <722  in  the  right- 
side  one  cam  be  observed,  while  the  corresponding  strains  61 1  and  £22  are  linearly 
vaurying.  Tramsverse  strains  £22  or  £ii  aure  present,  while  stresses  <r22  or  <7xx  ame  zero. 
The  transverse  shear  stress  <ri2  and  straun  £12  are  agaun  zero. 

Even  if,  in  this  exaimple,  the  rectilineau’  shape  of  the  beam  walls  prevents  any 
interlaminair  stress  to  au'ise,  it  is  expected  that  the  presence  of  a  transverse  stress 
state  would  induce  interlaminair  stresses  in  the  caae  of  curved  walls.  Indeed,  inter¬ 
laminar  stresses  au'e  adso  present  in  this  example  in  the  comer  areas,  but  a  more 
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Figure  22:  Box-beam,  upp>er-side  wall.  Strain  under  axial  load. 
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le.e  ic.e  a.e 
flb/inS) 


Strosa  Component  S23  (lb, 


Figure  23:  Box-be2an,  upper-side  wall.  Stress  under  axial  load. 
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Figure  26;  Box-beam,  upper-side  wall.  Strmn  under  torque. 
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Figure  27:  Box- beam,  upper-side  wall.  Stress  under  torque. 
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Figure  30:  The  curved  cross-ply  sp>ecimen. 


Geometric  dimensions 

Elastic  properties 

R  =  4.5  mm 
h  =  1.0  mm 
b  =  6.0  mm 

stacking  sequence:  [0/90], 

Eu  =  1407  Akg /mm} 

E22  =  E33  =  lA7Skg/mm} 

^13  =  =  '21  1/23  =  -3 

Gi2  =  Gi3  =  G23  —  598kg/mm^ 

Table  8:  Cross-ply  specimen.  Geometric  and  material  data. 


detailed  mesh  would  be  necessarj'  in  order  to  correctly  capture  them,  and  that  was 
beyond  the  scope  of  this  test  case. 

5.3  A  Curved  Cross-ply  Specimen 

A  comparison  with  a  fully  three-dimensional  approach  has  been  attempted  ana¬ 
lyzing  a  curved  composite  specimen  made  of  four  symmetric  graphite-epoxy  cross- 
ply  layers.  This  specimen,  see  Fig.  30  and  Table  8,  is  characterized  by  a  strong 
curvature  (i2//i=4.5,  il/6=. 75)  and  extends  over  180°:  the  mid  cross-section  is  con¬ 
sidered  when  the  specimen  is  loaded  as  shown  with  a  resultemt  force  of  1000  kg. 

The  three-dimensional  mesh  is  made  with  20-node  isoparametric  bride  elements 
to  model  each  ply;  because  of  symmetry,  a  quarter  of  the  sp>ecimen  has  been  modeled 
for  a  total  of  432  elements,  see  Fig.  31.  For  the  two-dimenaonal  analysis  of  the 
cross-section  three  meshes  have  been  modeled  2is  shown  in  Fig.  32:  a  coarse  one 
equivalent  to  the  3-D  model  as  fzu’  as  the  section  is  concerned  (4x6  parabolic  plane 
elements  for  half  a  section);  a  medium  one  with  2  elements  for  each  layer  and  still  6 
elements  on  the  section  half-span;  and  a  refined  one  with  3  elements  for  each  layer 
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Figtire  32:  Cross-ply  specimen.  2-D  meshes  of  a  half  of  the  cross-section:  (a)  coarse, 
(6)  medium,  (c)  fine. 


61 


Mesh 

Nodes 

DOFs 

Elements 

3-D 

2397 

7191 

432 

2-D  coarse 

93 

285 

24 

2-D  medium 

159 

483 

96 

2-D  fine 

481 

1449 

144 

Table  9:  Cross-ply  specimen.  Sizes  of  the  finite  element  models. 


and  12  elements  on  the  section  half-span.  Using  these  three  models  ensures  reaching 
convergence  of  restilts.  The  sizes  of  the  finite  element  models  are  given  in  Table  9: 
these  data  stress  once  more  the  effectiveness  of  the  present  approach. 

Results  are  compared  in  terms  of  stress  and  strain  across  the  thickness  at  the 
center  of  the  specimen.  The  stress  state  near  the  end  sides  of  the  specimen  section 
is  under  free  edge  effect  and  it  is  not  reported  here:  interlaminar  layers  should  have 
been  modeled  in  order  to  study  these  effects  and  this  was  not  the  purpose  of  this 
example.  As  shown  in  Fig.  33,  the  three  2-D  models  yield  the  same  results,  showing 
that  even  the  coarse  mesh  is  adequate.  Moreover,  both  the  2-D  and  3-D  analyses  are 
able  to  handle  the  discontinuities  in  the  material  elastic  properties:  the  stress  <733 
and  strain  €33  normal  to  the  cross-section  are  in  good  agreement.  On  the  contrary, 
the  two  transverse  stress  and  strain  components  differ  in  magnitude  in  2-D  and  3-D 
models.  This  discrepancy  would  probably  reduce  with  a  more  refined  discretization 
along  the  curved  beam  in  the  3-D  model:  in  fact,  the  displacement  representation 
along  the  beam  axis  in  the  present  approach  must  be  considered  exact 

5.4  A  Curved  Composite  I-Beam 

The  last  numerical  example  refers  to  a  curved  I-shaped  beam:  this  beam  was  first 
Euialyzed,  using  a  different  approach,  by  Peck  (1991)  and  it  constitutes  a  good 
numerical  test  for  the  present  method.  The  analysis  target  is  to  calculate  the  axial 
force  and  the  curling  moment  acting  up>on  the  flange  elements  of  the  beam  loaded 
independently  with  a  bending  moment  and  an  axizd  load.  A  peculiarity  of  this 
example  is  the  high  shear  force  transferred  between  the  flanges  and  web  due  to 
the  curvature,  even  under  bending  moment  zmd  axial  loads:  therefore  the  curling 
moment  of  the  flanges  becomes  a  significant  test  qucmtity. 

The  geometric  data  of  this  example  are  reported  in  Fig.  34,  and  sign  conventions 
in  Fig.  35.  Four  different  laminations  have  been  tested,  as  described  in  Table  10:  the 
elastic  properties  of  the  material,  graphite-epoxy  GR/EP  T300/5208,  are  reported 
in  Table  11.  The  finite  element  discretization  of  the  cross-section  is  stimmarized 
as  follows,  see  Fig.  36:  8x2  elements  across  the  thickness  of  the  web  (two  plane 
elements  averaging  each  one  6  layers),  and  12  x  22  elements  across  the  thickness  of 
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Coor^Mte  acroM  the  thielmaw  (ia)  Coer4iMte  acroaa  the  thiekaeae  (in) 


Figure  33:  Cross-ply  specimen.  Some  stress  and  strain  components  across  the 
thickness  at  the  center  of  the  specimen. 
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4 


Figure  34:  The  curved  composite  I-shaped  beam  (dimensions  are  in  inches). 


Nil  Nil 


Figure  35:  I-beam.  Conventions  for  pKJsitive  internal  forces. 


Table  10:  I-beam.  Stacking  sequences. 
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Table  11:  I-beam.  Elastic  properties. 


Figure  36:  I-beam.  Mesh  of  the  cross-section. 


65 


Case 

Axial  rigidity 

(lb) 

Bending  rigidity 
(lb  in.-2) 

1 

5.2815  xl0« 

1.7305  xlO^ 

2 

5.8488  xl0« 

1.9591  xlO^ 

3 

5.4504  xlO* 

1.8466  xlO^ 

4 

1.5726  xlO® 

3.8004  xlO® 

Table  12:  I-beam.  Cross-sectional  rigidities. 


each  flange  (a  plane  element  for  each  layer).  All  the  elements  are  8-node  isoparamet¬ 
ric  elements.  The  total  number  of  degrees  of  freedom  is  5343:  a  three-dimensional 
discretization  of  the  whole  beam  giving  such  detail  would  be  computationzJly  ex¬ 
pensive. 

Axial  and  bending  rigidities  in  the  analyzed  cases  are  reported  in  Table  12:  it  has 
been  observed  that  in  the  first  three  cases,  zdthough  the  stiffnesses  axe  comparable 
in  magnitude,  the  distribution  of  the  axial  force  and  the  curling  moment  on  the 
flanges  are  quite  different. 

The  deformed  cross-sections,  the  axial  stress  and  the  curling  moment  in  the 
flanges  due  to  tension  and  bending  moment  for  Ceise  1  are  shown  in  Figs.  37,  38 
and  39.  The  axial  tension  is  far  from  being  constant  along  the  flange  span,  showing 
that  this  kind  of  cross-section  is  not  an  optimum  with  respect  to  the  maximum  stress 
design  criterium.  This  is  mainly  due  to  the  interaction  of  the  curling  moment  with 
the  axial  stress.  Some  strain  components  are  plotted  in  Fig.  40  across  the  upper 
flange  thickness  (coordinates  X  and  Y  are  measured  with  reference  to  the  frame 
depicted  in  Fig.  35).  It  can  be  observed  that  while  the  extensional  strain  appeeu-s 
to  be  linear,  the  in-plane  shear  strain  does  not  behave  in  such  a  way.  This  fact 
suggests  that,  for  composite  sections,  it  czm  be  wrong  to  adopt  a  linear  variation  of 
the  displacement  through  the  thickness,  which  is  usual  in  simpler  approaches. 

In  Case  2  the  flanges  are  made  of  uniaxial  laminate  with  fibers  eilong  the  beam 
axis,  thus  increasing  the  section  axial  zmd  bending  stiffness.  Axial  forces  in  the 
flanges  come  out  to  be  about  twice  than  Case  1  (see  Figs.  42  and  43),  but  they  can 
not  spread  along  the  flange  width,  due  to  the  flange’s  poor  stiffness  to  withstand 
the  curling  moment,  which  is  lower  than  Case  1.  That  becomes  evident  from  the 
deformed  shape  shown  in  Fig.  41  and  is  confirmed  by  the  strains  plotted  in  Fig.  44: 
the  flanges  do  not  deform  significantly  near  their  free  edge  and  thus  do  not  cooperate 
to  the  structural  strength.  This  analysis  stresses  that  the  leunination  of  Case  2  is 
not  a  good  choice  as  it  would  be  for  a  straight  beam. 

Case  3  is  similar  to  Case  2  but  four  layers  have  fibers  along  the  flange  span 
thus  giving  some  capability  to  withstand  the  curling  moment.  That  yields  a  more 
uniform  distribution  along  the  flange  span  of  the  axial  stress,  see  Figs.  46  and  47. 
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a  b 


Figure  37:  I-beam,  Case  1.  Deformed  cross-section  under  (a)  tension,  (b)  bending 
moment. 

Nevertheless,  the  improvement  is  poor  with  respect  to  Ceise  2,  since  the  transverse 
plies  are  located  at  the  middle  of  the  flange  thickness:  it  can  be  seen  from  the 
deformation  results  in  Figs.  45  and  48  that  a  wide  area  of  the  fleinge  does  not  work 
significantly.  After  all,  the  axied  and  bending  stiffness  of  the  section  is  reduced  with 
respect  to  Case  2. 

In  Case  4  all  the  plies  lie  along  the  flange  span:  this  is  the  best  configuration 
as  regards  to  withstamd  the  curling  moment,  but  it  lowers  the  axial  strength  of  the 
flanges.  As  a  consequence,  the  flanges  do  not  quite  work  at  all  under  tension  and 
bending  and  the  web  is  called  to  carry  axial  and  bending  loads.  The  results  for 
Case  4  are  shown  in  Figs.  49  to  52. 

The  present  results  were  compared  with  the  three-dimensional  analysis  by  Peck 
(1991)  on  the  beam  model  shown  in  Fig.  53.  The  MSC-NASTRAN  mesh  consists  of 
350  8-node  plate  elements  and  eunounts  to  9000  DOFs.  In  the  central  cross-section 
the  applied  losids  give  a  compression  axial  force  of  1593  lb  and  a  bending  moment 
of  9558  lb  in.  The  axial  stress  and  the  curling  moment  in  the  flanges  axe  shown  in 
Fig.  54  together  with  the  results  of  the  present  approach. 
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Figure  40:  I-beam,  Case  1.  Some  strmn  components  across  the  thickness  of  the 
upper  flange  due  to  tension. 
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Figure  44:  I-beam,  Case  2.  Some  strain  components  across  the  thickness  of  the 
upper  flange  due  to  tension. 
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Figure  45:  I-beam,  Case  3.  Deformed  cross-section  under  (a)  tension,  (6)  bending 
moment. 
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Figure  48:  I-beam,  Case  3.  Some  strain  components  across  the  thickness  of  the 
upper  flange  due  to  tension. 
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Figure  49:  I-beam,  Case  4.  Deformed  cross-section  under  (a)  tension,  (6)  bending 
moment. 
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3-D  ANALYSIS 
PRESENT  APPROACH 


6  Conclusions 


The  method  developed  to  obtain  the  cross-sectional  constitutive  equations  of  natu¬ 
rally  curved  and  twisted  anisotropic  beams  is  consistent  with  the  three-dimensional 
theory  of  elasticity.  No  sp>ecific  simplifications  have  been  introduced  except  con¬ 
stancy  of  section,  curvature  and  twist,  which  are  commonly  used  in  beam  theory. 
The  potential  of  the  method  is  implied  in  the  fact  that  it  requires  a  two-dimensional 
approximation,  thus  achieving  a  degree  of  detail  which  is  considered  unfeasible  with 
a  three-dimensional  approach.  The  finite  element  einaJysis  of  the  beam  cross-section 
constitutes  a  good  step  towards  the  understanding  the  complexity  of  composite 
beams  and  a  useful  mathematical  tool  for  the  tailoring  of  real  designs. 

The  scope  of  the  research,  contained  in  the  proposal,  has  been  reached.  Never¬ 
theless,  additional  research  work  in  this  direction  should  be  done;  in  particular,  it 
is  in  the  authors’  opinion  that  it  would  be  extremely  useful  and  profitable  to  choose 
as  benchmark  some  significant  test  case  taken  from  real  design;  these  should  be 
submitted  to  different  research  centers  and  helicopter  industries  in  order  to  assess 
the  state  of  the  art  in  the  computation  of  composite  blades. 
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